]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ITS/UPGRADE/AliITSUSimulationPix.cxx
fix in calling of gaussian spread function
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUSimulationPix.cxx
index b5a0619b7fe2f19f05e848d98a0fd283b44edac5..cb6cc41d1d4893f527a6938e8fc5050ee84535ec 100644 (file)
@@ -1,8 +1,3 @@
-/*
- questions to experts: why RemoveDeadPixels should be called before FrompListToDigits ? 
-
-*/
 
 /**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
 
-#include <Riostream.h>
 #include <TGeoGlobalMagField.h>
 #include <TH1.h>
+#include <TH2.h>
 #include <TString.h>
 #include "AliITSU.h"
 #include "AliITSUDigitPix.h"
 #include "AliITSUHit.h"
-#include "AliITSUModule.h"
+#include "AliITSUChip.h"
 #include "AliITSUSensMap.h"
 #include "AliITSUCalibrationPix.h"
 #include "AliITSUSegmentationPix.h"
 #include "AliMathBase.h"
 #include "AliITSUSimuParam.h"
 #include "AliITSUSDigit.h"
+#include "AliITSUParamList.h"
+
+using std::cout;
+using std::endl;
+using namespace TMath;
 
 ClassImp(AliITSUSimulationPix)
 ////////////////////////////////////////////////////////////////////////
@@ -54,27 +54,33 @@ ClassImp(AliITSUSimulationPix)
 //
 //  AliITSUSimulationPix is to do the simulation of pixels
 //
+//  2013 Feb: Added MonoPix response and nois calculation al la MIMOSA32 (levente.molnar@cern.ch)
+//
 ////////////////////////////////////////////////////////////////////////
 
 //______________________________________________________________________
 AliITSUSimulationPix::AliITSUSimulationPix()
 :  fTanLorAng(0)
-  ,fStrobe(kTRUE)
-  ,fStrobeLenght(4)
-  ,fStrobePhase(-12.5e-9)
+  ,fGlobalChargeScale(1.0)
+  ,fSpread2DHisto(0)
+  ,fSpreadFun(0)
+  ,fROTimeFun(0)
 {
    // Default constructor.
+  SetUniqueID(AliITSUGeomTGeo::kChipTypePix);
 }
 
 //______________________________________________________________________
 AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap* map)
   :AliITSUSimulation(sim,map)
   ,fTanLorAng(0)
-  ,fStrobe(kTRUE)
-  ,fStrobeLenght(4)
-  ,fStrobePhase(-12.5e-9)
+  ,fGlobalChargeScale(1.0)
+  ,fSpread2DHisto(0)
+  ,fSpreadFun(0)
+  ,fROTimeFun(0)
 {
   // standard constructor
+  SetUniqueID(AliITSUGeomTGeo::kChipTypePix);
   Init();
 }
 
@@ -82,9 +88,10 @@ AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap*
 AliITSUSimulationPix::AliITSUSimulationPix(const AliITSUSimulationPix &s) 
   :AliITSUSimulation(s)
   ,fTanLorAng(s.fTanLorAng)
-  ,fStrobe(s.fStrobe)
-  ,fStrobeLenght(s.fStrobeLenght)
-  ,fStrobePhase(s.fStrobePhase)
+  ,fGlobalChargeScale(s.fGlobalChargeScale)
+  ,fSpread2DHisto(s.fSpread2DHisto)
+  ,fSpreadFun(s.fSpreadFun)
+  ,fROTimeFun(s.fROTimeFun)
 {
   //     Copy Constructor
 }
@@ -103,9 +110,12 @@ AliITSUSimulationPix& AliITSUSimulationPix::operator=(const AliITSUSimulationPix
   //    Assignment operator
   if (&s == this) return *this;
   AliITSUSimulation::operator=(s);
-  fStrobe       = s.fStrobe;
-  fStrobeLenght = s.fStrobeLenght;
-  fStrobePhase  = s.fStrobePhase;
+  fSpread2DHisto = s.fSpread2DHisto;
+  //
+  fGlobalChargeScale = s.fGlobalChargeScale;
+  fSpreadFun    = s.fSpreadFun;
+  fROTimeFun    = s.fROTimeFun;
+  //
   return *this;
 }
 
@@ -138,29 +148,22 @@ Bool_t AliITSUSimulationPix::SetTanLorAngle(Double_t weightHole)
   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
   if (!fld) AliFatal("The field is not initialized");
   Double_t bz = fld->SolenoidField();
-  fTanLorAng = TMath::Tan(weightHole*fSimuParam->LorentzAngleHole(bz) +
+  fTanLorAng = Tan(weightHole*fSimuParam->LorentzAngleHole(bz) +
                          weightEle*fSimuParam->LorentzAngleElectron(bz));
   return kTRUE;
 }
 
 //_____________________________________________________________________
-void AliITSUSimulationPix::SDigitiseModule(AliITSUModule *mod, Int_t /*mask*/,Int_t event, AliITSsegmentation* seg)
+void AliITSUSimulationPix::SDigitiseChip()
 {
   //  This function begins the work of creating S-Digits.
-  if (!(mod->GetNHits())) {
-    AliDebug(1,Form("In event %d module %d there are %d hits returning.",
-                   event, mod->GetIndex(),mod->GetNHits()));
-    return;
-  } 
-  //
-  if (fStrobe) if (event != fEvent) GenerateStrobePhase(); 
-  InitSimulationModule(mod->GetIndex(), event, seg );
-  // Hits2SDigits(mod);
-  Hits2SDigitsFast(mod);
-  if (fSimuParam->GetPixAddNoisyFlag())   AddNoisyPixels();
-  if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();  
+    
+  AliDebug(10,Form("In event %d chip %d there are %d hits", fEvent, fChip->GetIndex(),fChip->GetNHits()));       
+  if (fChip->GetNHits()) Hits2SDigitsFast();
+  if (!fSensMap->GetEntries()) return;
   WriteSDigits();
   ClearMap();
+  //
 }
 
 //______________________________________________________________________
@@ -178,7 +181,7 @@ void AliITSUSimulationPix::WriteSDigits()
 }
 
 //______________________________________________________________________
-void AliITSUSimulationPix::FinishSDigitiseModule()
+void AliITSUSimulationPix::FinishSDigitiseChip()
 {
    //  This function calls SDigitsToDigits which creates Digits from SDigits
   FrompListToDigits();
@@ -187,80 +190,61 @@ void AliITSUSimulationPix::FinishSDigitiseModule()
 }
 
 //______________________________________________________________________
-void AliITSUSimulationPix::DigitiseModule(AliITSUModule *mod,Int_t /*mask*/, Int_t event, AliITSsegmentation* seg)
+void AliITSUSimulationPix::DigitiseChip()
 {
   //  This function creates Digits straight from the hits and then adds
   //  electronic noise to the digits before adding them to pList
   //  Each of the input variables is passed along to Hits2SDigits
   //
-  if (fStrobe) if (event != fEvent) GenerateStrobePhase();
-  InitSimulationModule( mod->GetIndex(), event, seg );
-  // Hits2SDigits(mod);
-  Hits2SDigitsFast(mod);
-  //
-  if (fSimuParam->GetPixAddNoisyFlag())   AddNoisyPixels();
-  if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();
-  FrompListToDigits();
-  ClearMap();
+  // pick charge spread function
+  Hits2SDigitsFast();
+  FinishSDigitiseChip();
 }
 
 //______________________________________________________________________
-void AliITSUSimulationPix::Hits2SDigits(AliITSUModule *mod)
+void AliITSUSimulationPix::Hits2SDigits()
 {
   // Does the charge distributions using Gaussian diffusion charge charing.
-  const Double_t kmictocm = 1.0e-4; // convert microns to cm.
-  const Double_t kcmtomic = 1.e4;
-  const Double_t kBunchLenght = 25e-9; // LHC clock
-  Int_t nhits = mod->GetNHits();
+  Int_t nhits = fChip->GetNHits();
   if (!nhits) return;
   //
   Int_t h,ix,iz,i;
   Int_t idtrack;
   Float_t x,y,z;  // keep coordinates float (required by AliSegmentation)
-  Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0,ld=0.0;
-  Double_t t,tp,st,dt=0.2,el,sig,sigx,sigz,fda;
-  Double_t thick = 0.5*kmictocm*fSeg->Dy();  // Half Thickness
-  fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
+  Double_t tof,x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0;
+  Double_t t,tp,st,dt=0.2,el;
+  Double_t thick = 0.5*fSeg->Dy();  // Half Thickness
+
   //
   for (h=0;h<nhits;h++) {
-    if (fStrobe && 
-       ((mod->GetHit(h)->GetTOF()<fStrobePhase) || 
-        (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
-       ) continue;
     //
-    if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
-    st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
+    if (!fChip->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
+    st = Sqrt(x1*x1+y1*y1+z1*z1);
     if (st>0.0) {
-      st = (Double_t)((Int_t)(st*kcmtomic)); // number of microns
+      st = (Double_t)((Int_t)(st*1e4)); // number of microns
       if (st<=1.0) st = 1.0;
-      dt = 1.0/st;
+      dt = 1.0/st;               // RS TODO: do we need 1 micron steps?
+      double dy = dt*thick;
+      y = -0.5*dy;
       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;
+       y  += dy;
        z   = z0+z1*tp;
        if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
-       el  = dt * de / fSimuParam->GetGeVToCharge();
-       //
-       if (GetDebug(1)) if (el<=0.0) cout<<"el="<<el<<" dt="<<dt<<" de="<<de<<endl; // end if GetDebug
+       el  = fGlobalChargeScale * dt * de / fSimuParam->GetGeVToCharge();
        //
-       sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y)); 
-       sigx=sig;
-       sigz=sig*fda;
-       if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
-       SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
+       if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
+       SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
       } // end for t
     } else { // st == 0.0 deposit it at this point
       x   = x0;
-      y   = y0;
+      y   = y0 + 0.5*thick;
       z   = z0;
       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
-      el  = de / fSimuParam->GetGeVToCharge();
-      sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
-      sigx = sig;
-      sigz = sig*fda;
-      if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
-      SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
+      el  = fGlobalChargeScale * de / fSimuParam->GetGeVToCharge();
+      if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
+      SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
     } // end if st>0.0    
   } // Loop over all hits h
   //
@@ -268,18 +252,20 @@ void AliITSUSimulationPix::Hits2SDigits(AliITSUModule *mod)
   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
   AliITSUSDigit* dg = 0;
   switch (fSimuParam->GetPixCouplingOption()) {
+  case AliITSUSimuParam::kNoCouplingPix : 
+    break;
   case AliITSUSimuParam::kNewCouplingPix :
     for (i=nd;i--;) {
       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
       if (fSensMap->IsDisabled(dg)) continue;
-      SetCoupling(dg,idtrack,h);
+      SetCoupling(dg);
     } 
     break;
   case AliITSUSimuParam::kOldCouplingPix:
     for (i=nd;i--;) {
       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
       if (fSensMap->IsDisabled(dg)) continue;
-      SetCouplingOld(dg,idtrack,h);
+      SetCouplingOld(dg);
     } 
     break;
   default:
@@ -290,72 +276,56 @@ void AliITSUSimulationPix::Hits2SDigits(AliITSUModule *mod)
 }
 
 //______________________________________________________________________
-void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod)
+void AliITSUSimulationPix::Hits2SDigitsFast()
 {
   // Does the charge distributions using Gaussian diffusion charge charing.    // Inputs:
-  //    AliITSUModule *mod  Pointer to this module
-  // Output:
-  //    none.
-  // Return:
-  //    none.
-  const Double_t kmictocm = 1.0e-4; // convert microns to cm.
-  const Int_t kn10=10;
-  const Double_t kti[kn10]={7.443716945e-3,2.166976971e-1,3.397047841e-1,
-                           4.325316833e-1,4.869532643e-1,5.130467358e-1,
-                           5.674683167e-1,6.602952159e-1,7.833023029e-1,
-                           9.255628306e-1};
-  const Double_t kwi[kn10]={1.477621124e-1,1.346333597e-1,1.095431813e-1,
-                           7.472567455e-2,3.333567215e-2,3.333567215e-2,
-                           7.472567455e-2,1.095431813e-1,1.346333597e-1,
-                           1.477621124e-1};
-  const Double_t kBunchLenght = 25e-9; // LHC clock
-  TObjArray *hits = mod->GetHits();
+  //    AliITSUChip *mod  Pointer to this chip
+  //
+  TObjArray *hits = fChip->GetHits();
   Int_t nhits = hits->GetEntriesFast();
   if (nhits<=0) return;
   //
   Int_t h,ix,iz,i;
   Int_t idtrack;
   Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
-  Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0; 
-  Double_t t,st,el,sig,sigx,sigz,fda,de=0.0,ld=0.0;
-  Double_t thick = 0.5*kmictocm*fSeg->Dy();  // Half thickness
-  fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
+  Double_t tof,x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0; 
+  Double_t step,st,el,de=0.0;
+  Double_t minDim = Min(fSeg->Dpx(1),fSeg->Dpz(1)); // RStmp: smallest pitch
+  Double_t thick = fSeg->Dy();
   //
   for (h=0;h<nhits;h++) {
-    // Check if the hit is inside readout window
-    if (fStrobe && 
-       ((mod->GetHit(h)->GetTOF()<fStrobePhase) || 
-        (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
-       ) continue;
     //
-    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) 
-      for (i=0;i<kn10;i++) { // Integrate over t
-       t   = kti[i];
-       x   = x0+x1*t;
-       y   = y0+y1*t;
-       z   = z0+z1*t;
+    if (!fChip->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
+    //
+    st = Sqrt(x1*x1+y1*y1+z1*z1); 
+    if (st>0.0) {
+      int np = int(1.5*st/minDim);  //RStmp: inject the points in such a way that there is ~1.5 point per cell
+      np = TMath::Max(1.0*np,fResponseParam->GetParameter(kSpreadFunMinSteps));   
+      AliDebug(10,Form(" Number of charge injection steps is set to %d ",np));   
+      double dstep = 1./np;
+      double dy = dstep*thick;
+      y = -0.5*dy;
+      step = -0.5*dstep;
+      for (i=0;i<np;i++) {          //RStmp Integrate over t
+       //      for (i=0;i<kn10;i++) { // Integrate over t
+       step  += dstep;  // RStmp kti[i];
+       x   = x0+x1*step;
+       y  += dy;
+       z   = z0+z1*step;
        if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
-       el  = kwi[i]*de/fSimuParam->GetGeVToCharge();
-       sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
-       sigx=sig;
-       sigz=sig*fda;
-       if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
-       SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
-       //                cout << "sigx sigz " << sigx << " " << sigz << endl; // dom
+       el  = fGlobalChargeScale * dstep * de/fSimuParam->GetGeVToCharge();
+       if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
+       SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
       } // end for i // End Integrate over t
+    }
     else { // st == 0.0 deposit it at this point
       x   = x0;
-      y   = y0;
+      y   = y0+0.5*thick;
       z   = z0;
       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
-      el  = de / fSimuParam->GetGeVToCharge();
-      sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
-      sigx=sig;
-      sigz=sig*fda;
-      if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
-      SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
+      el  = fGlobalChargeScale * de / fSimuParam->GetGeVToCharge();
+      if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
+      SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
     } // end if st>0.0
     
   } // Loop over all hits h
@@ -364,17 +334,19 @@ void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod)
   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
   AliITSUSDigit* dg = 0;
   switch (fSimuParam->GetPixCouplingOption()) {
+  case AliITSUSimuParam::kNoCouplingPix : 
+    break;
   case AliITSUSimuParam::kNewCouplingPix :
     for (i=nd;i--;) {
       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
       if (fSensMap->IsDisabled(dg)) continue;
-      SetCoupling(dg,idtrack,h);
+      SetCoupling(dg);
     } 
   case AliITSUSimuParam::kOldCouplingPix:
     for (i=nd;i--;) {
       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
       if (fSensMap->IsDisabled(dg)) continue;
-      SetCouplingOld(dg,idtrack,h);
+      SetCouplingOld(dg);
     } 
     break;
   default:
@@ -384,178 +356,166 @@ void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod)
 }
 
 //______________________________________________________________________
-void AliITSUSimulationPix::SpreadCharge(Double_t x0,Double_t z0,
-                                         Int_t ix0,Int_t iz0,
-                                         Double_t el,Double_t sig,Double_t ld,
-                                         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)
-   // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
-   // 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 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)
-   //    Double_t ld   lorentz drift in x 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;  // keep coordinates float (required by AliSegmentation)
-   Double_t s,sp,x1,x2,z1,z2; 
-   //
-   if (GetDebug(2)) AliInfo(Form("(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.
-     UpdateMapSignal(iz0,ix0,t,hi,el);
-     return;
-   } // end if
-   sp = 1.0/(sig*kRoot2);
-   ixs = TMath::Max(-knx+ix0,0);
-   ixe = TMath::Min(knx+ix0,fSeg->Npx()-1);
-   izs = TMath::Max(-knz+iz0,0);
-   ize = TMath::Min(knz+iz0,fSeg->Npz()-1);
-   for (ix=ixs;ix<=ixe;ix++) 
-     for (iz=izs;iz<=ize;iz++) {
-       fSeg->DetToLocal(ix,iz,x,z); // pixel center
-       double dxi = 0.5*kmictocm*fSeg->Dpx(ix);
-       double dzi = 0.5*kmictocm*fSeg->Dpz(iz);
-       x1  = x;
-       z1  = z;
-       x2  = x1 + dxi; // Upper
-       x1 -= dxi;  // Lower
-       z2  = z1 + dzi; // Upper
-       z1 -= dzi;  // Lower
-       x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
-       x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
-       z1 -= z0; // Distance from where track traveled
-       z2 -= z0; // Distance from where track traveled
-       s   = el*0.25; // Correction based on definision of Erfc
-       s  *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(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  *= AliMathBase::ErfcFast(sp*z1) - AliMathBase::ErfcFast(sp*z2);
-       if (GetDebug(3)) cout<<" sp*z1="<<sp*z1<<" sp*z2="<<sp*z2<<" s="<<s<< endl; // end if GetDebug
-       if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
-     } // end for ix, iz
-   //
-}
-
-//______________________________________________________________________
-void AliITSUSimulationPix::SpreadChargeAsym(Double_t x0,Double_t z0,
-                                           Int_t ix0,Int_t iz0,
-                                           Double_t el,Double_t sigx,Double_t sigz,
-                                           Double_t ld,Int_t t,Int_t hi)
+void AliITSUSimulationPix::SpreadCharge2D(Double_t x0,Double_t z0, Double_t dy, Int_t ix0,Int_t iz0,
+                                         Double_t el, Double_t tof, Int_t tID, Int_t hID)
 {
   // Spreads the charge over neighboring cells. Assume charge is distributed
   // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
   // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
-  // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
   // 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 z0   z position of point where charge is liberated
+  //    Double_t x0   x position of point where charge is liberated (local)
+  //    Double_t z0   z position of point where charge is liberated (local)
+  //    Double_t dy   distance from the entrance surface (diffusion sigma may depend on it)
   //    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 sigx Sigma difusion along x for this step (y0 dependent)
   //    Double_t sigz Sigma difusion along z for this step (z0 dependent)
-  //    Double_t ld   lorentz drift in x for this stip (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 = 3; // RS: TO TUNE
-  const Double_t kRoot2 = 1.414213562; // Sqrt(2).
-  const Double_t kmictocm = 1.0e-4; // convert microns to cm.
+  //    Int_t    tID  track number
+  //    Int_t    hID  hit "hit" index number
+  //
   Int_t ix,iz,ixs,ixe,izs,ize;
   Float_t x,z;   // keep coordinates float (required by AliSegmentation)
-  Double_t s,spx,spz,x1,x2,z1,z2; 
-  //
-  if (GetDebug(2)) AliInfo(Form("SpreadChargeAsym","(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,sigx=%e, sigz=%e, t=%d,i=%d,ld=%e)",x0,z0,ix0,iz0,el,sigx,sigz,t,hi,ld));
-  if (sigx<=0.0 || sigz<=0.0) { // if sig<=0 No diffusion to simulate.
-    UpdateMapSignal(iz0,ix0,t,hi,el);
-    if (GetDebug(2)) {
-      cout << "sigx<=0.0=" << sigx << endl;
-      cout << "sigz<=0.0=" << sigz << endl;
-    } // end if GetDebug
-    return;
-  } // end if
-  spx = 1.0/(sigx*kRoot2);     
-  spz = 1.0/(sigz*kRoot2);
-  if (GetDebug(2)) {
-    cout << "sigx=" << sigx << " spx=" << spx << endl;
-    cout << "sigz=" << sigz << " spz=" << spz << endl;
-  } // end if GetDebug
-  ixs = TMath::Max(-knx+ix0,0);
-  ixe = TMath::Min(knx+ix0,fSeg->Npx()-1);
-  izs = TMath::Max(-knz+iz0,0);
-  ize = TMath::Min(knz+iz0,fSeg->Npz()-1);
+  Float_t xdioshift = 0 , zdioshift = 0 ;
+  Double_t s,dtIn[kNDtSpread]; // data transfered to spread function for integral calculation
+  //
+  if (GetDebug(2)) AliInfo(Form("(x0=%e,z0=%e,dy=%e, ix0=%d,iz0=%d,el=%e,tID=%d,hID=%d)",x0,z0,dy,ix0,iz0,el,tID,hID));
+  //
+  Double_t &x1 = dtIn[kCellX1];
+  Double_t &x2 = dtIn[kCellX2];
+  Double_t &z1 = dtIn[kCellZ1];
+  Double_t &z2 = dtIn[kCellZ2];
+  //
+  int nx = GetResponseParam()->GetParameter(kSpreadFunParamNXoffs);
+  int nz = GetResponseParam()->GetParameter(kSpreadFunParamNZoffs);
+  //
+  dtIn[kCellYDepth]  = dy;
+  ixs = Max(-nx+ix0,0);
+  ixe = Min( nx+ix0,fSeg->Npx()-1);
+  izs = Max(-nz+iz0,0);
+  ize = Min( nz+iz0,fSeg->Npz()-1);
   for (ix=ixs;ix<=ixe;ix++) 
     for (iz=izs;iz<=ize;iz++) {
+      //
+      // Check if the hit is inside readout window
+      int cycleRO = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof);
+      if (Abs(cycleRO)>kMaxROCycleAccept) continue;
+      //
       fSeg->DetToLocal(ix,iz,x,z); // pixel center
-      double dxi = 0.5*kmictocm*fSeg->Dpx(ix);
-      double dzi = 0.5*kmictocm*fSeg->Dpz(iz);
-      x1  = x;
-      z1  = z;
+      xdioshift = zdioshift = 0;
+      double dxi = fSeg->Dpx(ix);
+      double dzi = fSeg->Dpz(iz);
+      CalcDiodeShiftInPixel(ix,iz,xdioshift,zdioshift);    // Check and apply diode shift if needed
+      xdioshift *= dxi;
+      zdioshift *= dzi;
+      dxi *= 0.5;
+      dzi *= 0.5;
+      //      printf("DShift: %d %d -> %.4f %.4f\n",ix,iz,xdioshift,zdioshift);
+      x1  = (x + xdioshift) - x0;   // calculate distance of cell boundaries from injection center
+      z1  = (z + zdioshift) - z0;
       x2  = x1 + dxi; // Upper
       x1 -= dxi;      // Lower
       z2  = z1 + dzi; // Upper
       z1 -= dzi;      // Lower
-      x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
-      x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
-      z1 -= z0; // Distance from where track traveled
-      z2 -= z0; // Distance from where track traveled
-      s   = el*0.25; // Correction based on definision of Erfc
-      s  *= AliMathBase::ErfcFast(spx*x1) - AliMathBase::ErfcFast(spx*x2);
-      if (GetDebug(3)) cout <<"el="<<el<<" ix0="<<ix0<<" ix="<<ix<<" x0="<<x<<" iz0="<<iz0<<" iz="<<iz<<" z0="<<z<<" spx*x1="<<spx*x1<<" spx*x2="<<spx*x2<<" s="<<s; // end if GetDebug
-      s  *= AliMathBase::ErfcFast(spz*z1) - AliMathBase::ErfcFast(spz*z2);
-      if (GetDebug(3)) cout<<" spz*z1="<<spz*z1<<" spz*z2="<<spz*z2<<" s="<<s<< endl; // end if GetDebug
-      if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
+      s   = el* (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fSpreadFun)(dtIn);
+      if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,tID,hID,s,cycleRO);
     } // end for ix, iz
   //
 }
 
+//______________________________________________________________________
+Double_t AliITSUSimulationPix::SpreadFunDoubleGauss2D(const Double_t *dtIn)
+{
+  // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2] 
+  // and Z=dtIn[kCellZ1]:dtIn[kCellZ2] 
+  // The spread function is assumed to be double gaussian in 2D
+  // Parameters should be: mean0,sigma0, mean1,sigma1, relative strenght of 2nd gaussian wrt 1st one
+  //
+  // 1st gaussian
+  double intg1 = GausInt2D(fResponseParam->GetParameter(kG2SigX0),  // sigX
+                          dtIn[kCellX1]-fResponseParam->GetParameter(kG2MeanX0),      // x1-xmean
+                          dtIn[kCellX2]-fResponseParam->GetParameter(kG2MeanX0),      // x2-xmean
+                          fResponseParam->GetParameter(kG2SigZ0),  // sigZ
+                          dtIn[kCellZ1]-fResponseParam->GetParameter(kG2MeanZ0),    // z1-zmean
+                          dtIn[kCellZ2]-fResponseParam->GetParameter(kG2MeanZ0));   // z2-zmean
+  // 2nd gaussian
+  double intg2 = GausInt2D(fResponseParam->GetParameter(kG2SigX1),  // sigX
+                          dtIn[kCellX1]-fResponseParam->GetParameter(kG2MeanX1),      // x1-xmean
+                          dtIn[kCellX2]-fResponseParam->GetParameter(kG2MeanX1),      // x2-xmean
+                          fResponseParam->GetParameter(kG2SigZ1),  // sigZ
+                          dtIn[kCellZ1]-fResponseParam->GetParameter(kG2MeanZ1),    // z1-zmean
+                          dtIn[kCellZ2]-fResponseParam->GetParameter(kG2MeanZ1));   // z2-zmean
+  double scl = fResponseParam->GetParameter(kG2ScaleG2);
+  return (intg1+intg2*scl)/(1+scl);
+  //
+} 
+
+
+//______________________________________________________________________
+Double_t AliITSUSimulationPix::SpreadFrom2DHisto(const Double_t *dtIn)
+{
+  // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2] 
+  // and Z=dtIn[kCellZ1]:dtIn[kCellZ2] 
+  // The spread function integral is taken from fSpread2DHisto extracted from the sensor response parameters
+  // list in the method SetResponseParam. The histo must return the fraction of charge integrates in the 
+  // cell whose center is on the distance X=(dtIn[kCellX1]+dtIn[kCellX2])/2 and Z=(dtIn[kCellZ1]+dtIn[kCellZ2])/2
+  // from the injection point.
+  //
+  Double_t qpixfrac = 0;  
+  Double_t xintp = 1e4*(dtIn[kCellX1]+dtIn[kCellX2])/2.0;
+  Double_t zintp = 1e4*(dtIn[kCellZ1]+dtIn[kCellZ2])/2.0;
+  //    
+  qpixfrac =  fSpread2DHisto->Interpolate(xintp,zintp); //the PSF map is given in um but the dtIn is in cm so we need to convert it 
+  //
+  return qpixfrac;  
+}
+
+//______________________________________________________________________
+Double_t AliITSUSimulationPix::SpreadFunGauss2D(const Double_t *dtIn)
+{
+  // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2] 
+  // and Z=dtIn[kCellZ1]:dtIn[kCellZ2] 
+  // The spread function is assumed to be gaussian in 2D
+  // Parameters should be: mean0,sigma0
+  return GausInt2D(fResponseParam->GetParameter(kG1SigX),  // sigX
+                  dtIn[kCellX1]-fResponseParam->GetParameter(kG1MeanX),    // x1-xmean
+                  dtIn[kCellX2]-fResponseParam->GetParameter(kG1MeanX),    // x2-xmean
+                  fResponseParam->GetParameter(kG1SigZ),  // sigZ
+                  dtIn[kCellZ1]-fResponseParam->GetParameter(kG1MeanZ),    // z1-zmean
+                  dtIn[kCellZ2]-fResponseParam->GetParameter(kG1MeanZ));   // z2-zmean
+  //
+} 
+
 //______________________________________________________________________
 void AliITSUSimulationPix::RemoveDeadPixels() 
 {
-  // Removes dead pixels on each module (ladder)
-  // This should be called before going from sdigits to digits (FrompListToDigits)
+  // Removes dead pixels on each chip (ladder)
+  // This should be called before going from sdigits to digits (i.e. from FrompListToDigits)
   
   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibDead();
   if (!calObj) return;
   //
-  if (calObj->IsBad()) {ClearMap(); return;} // whole module is masked
+  if (calObj->IsBad()) {ClearMap(); return;} // whole chip is masked
   //
+  // prepare the list of r/o cycles seen
+  Char_t cyclesSeen[2*kMaxROCycleAccept+1];
+  int ncycles = 0;
+  for (int i=(2*kMaxROCycleAccept+1);i--;) if (fCyclesID[i]) cyclesSeen[ncycles++]=i-kMaxROCycleAccept;
+  
   // remove single bad pixels one by one
   int nsingle = calObj->GetNrBadSingle();
   UInt_t col,row;
+  Int_t cycle;
   for (int i=nsingle;i--;) {
     calObj->GetBadPixelSingle(i,row,col);
-    fSensMap->DeleteItem(col,row);
+    for (int icl=ncycles;icl--;) fSensMap->DeleteItem(col,row,cyclesSeen[icl]);
   }
   int nsd = fSensMap->GetEntriesUnsorted();
   for (int isd=nsd;isd--;) {
     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->AtUnsorted(isd);
     if (fSensMap->IsDisabled(sd)) continue;
-    fSensMap->GetMapIndex(sd->GetUniqueID(),col,row);
+    fSensMap->GetMapIndex(sd->GetUniqueID(),col,row,cycle);
     int chip = fSeg->GetChipFromChannel(0,col);
     //    if (calObj->IsChipMarkedBad(chip)) fSensMap->Disable(sd); // this will simple mark the hit as bad
     if (calObj->IsChipMarkedBad(chip)) fSensMap->DeleteItem(sd); // this will suppress hit in the sorted list
@@ -566,66 +526,62 @@ void AliITSUSimulationPix::RemoveDeadPixels()
 //______________________________________________________________________
 void AliITSUSimulationPix::AddNoisyPixels() 
 {
-  // Adds noisy pixels on each module (ladder)
-  // This should be called before going from sdigits to digits (FrompListToDigits)
+  // Adds noisy pixels on each chip (ladder)
+  // This should be called before going from sdigits to digits (i.e. FrompListToDigits)
   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibNoisy();
-  if (!calObj) return;
+  if (!calObj) { AliDebug(10,Form("  No Calib Object for Noise!!! ")); return;}
   for (Int_t i=calObj->GetNrBad(); i--;) UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i), 
-                                                       10*fSimuParam->GetPixThreshold(fModule));
+                                                       10*fSimuParam->GetPixThreshold(fChip->GetIndex()));
   //
 }
 
 //______________________________________________________________________
 void AliITSUSimulationPix::FrompListToDigits() 
 {
-  // add noise and electronics, perform the zero suppression and add the
-  // digit to the list
-  static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
-  UInt_t ix,iz;
-  Double_t sig;
-  const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
-  static AliITSUDigitPix dig;
-  // RS: in principle:
-  // 1) for every pixel w/o hit we have to generate a noise and activate the pixel if the noise exceeds the threshold. 
-  // 2) for every pixel with hit we should add random noise and check if the total signal exceeds the threshold
-  // With many channels this will be too time consuming, hence I do the following
-  // 1) Precalculate the probability that the nois alone will exceed the threshold. 
-  // 2) Chose randomly empty pixels according to this probability and apply the noise above threshold.
-  // 3) For pixels having a hits apply the usual noise and compare with threshold
+  // add noise and electronics, perform the zero suppression and add the digits to the list
   //
   // RS may use for ordered random sample generation dl.acm.org/ft_gateway.cfm?id=356313&type=pdf
   //
-  int maxInd = fSensMap->GetMaxIndex();
-  double minProb = 0.1/maxInd;
+  int nsd = fSensMap->GetEntriesUnsorted();  // sdigits added from the signal
   //
-  int nsd = fSensMap->GetEntries();
-  Int_t prevID=0,curID=0;
-  TArrayI ordSampleInd(100),ordSample(100);
+  // add different kinds of noise.
+  Bool_t addNoisy = fSimuParam->GetPixAddNoisyFlag() && (nsd>0 || fSimuParam->GetPixNoiseInAllMod()); // do we generate noise?
+  if (addNoisy) {
+    AddNoisyPixels();       // constantly noisy channels
+    AddRandomNoisePixels(0.0); // random noise: at the moment generate noise only for instance 0
+    nsd = fSensMap->GetEntriesUnsorted();
+  }
   //
-  double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(fModule);
-  fSimuParam->GetPixNoise(fModule, noiseSig, noiseMean);
-  probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold
+  if (nsd && fSimuParam->GetPixRemoveDeadFlag()) {
+    RemoveDeadPixels(); 
+    // note that here we shall use GetEntries instead of GetEntriesUnsorted since the 
+    // later operates on the array where the elements are not removed by flagged
+    nsd = fSensMap->GetEntries(); 
+  }
+  if (!nsd) return; // nothing to digitize
+  //
+  UInt_t row,col;
+  Int_t iCycle,modId = fChip->GetIndex();
+  Double_t sig;
+  const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
+  static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
+  static AliITSUDigitPix dig;
   //
   for (int i=0;i<nsd;i++) {
     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
     if (fSensMap->IsDisabled(sd)) continue;
-    curID = (int)sd->GetUniqueID();
-    //
-    if (probNoisy>minProb) { // generate randomly noisy pixels above the threshold, with ID's between previous hit and current
-      CreateNoisyDigits(prevID,curID,probNoisy, noiseSig, noiseMean);
-      prevID = curID+1;
-    }
     //
-    if ((sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(fModule)) continue;
-    if (TMath::Abs(sig)>2147483647.0) { //RS?
+    if ((sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(modId)) continue;
+    if (Abs(sig)>2147483647.0) { //RS?
       //PH 2147483647 is the max. integer
       //PH This apparently is a problem which needs investigation
       AliWarning(Form("Too big or too small signal value %f",sig));
     }
-    fSensMap->GetMapIndex(sd->GetUniqueID(),iz,ix);
-    dig.SetCoord1(iz);
-    dig.SetCoord2(ix);
-    dig.SetSignal(1);
+    fSensMap->GetMapIndex(sd->GetUniqueID(),col,row,iCycle);
+    dig.SetCoord1(col);
+    dig.SetCoord2(row);
+    dig.SetROCycle(iCycle);
+    dig.SetSignal((Int_t)sig);
     dig.SetSignalPix((Int_t)sig);
     int ntr = sd->GetNTracks();
     for (int j=0;j<ntr;j++) {
@@ -636,80 +592,60 @@ void AliITSUSimulationPix::FrompListToDigits()
       dig.SetTrack(j,-3);
       dig.SetHit(j,-1);
     }
-    aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix, &dig);
+    aliITS->AddSimDigit(AliITSUGeomTGeo::kChipTypePix, &dig);
   }
-  // if needed, add noisy pixels with id from last real hit to maxID
-  if (probNoisy>minProb) CreateNoisyDigits(prevID,maxInd,probNoisy, noiseSig, noiseMean);
   // 
 }
 
 //______________________________________________________________________
-Int_t AliITSUSimulationPix::CreateNoisyDigits(Int_t minID,Int_t maxID,double probNoisy, double noise, double base)
+Int_t AliITSUSimulationPix::AddRandomNoisePixels(Double_t tof) 
 {
-  // create random noisy digits above threshold within id range [minID,maxID[
-  // see FrompListToDigits for details
+  // create random noisy sdigits above threshold 
   //
-  static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
-  UInt_t ix,iz;
-  static AliITSUDigitPix dig;
+  int modId = fChip->GetIndex();
+  int npix = fSeg->GetNPads();
+  int ncand = gRandom->Poisson( npix*fSimuParam->GetPixFakeRate() );
+  if (ncand<1) return 0;
+  //
+  double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(modId);
+  fSimuParam->GetPixNoise(modId, noiseSig, noiseMean);
+  probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold
+  UInt_t row,col;
+  Int_t iCycle;
   static TArrayI ordSampleInd(100),ordSample(100); //RS!!! static is not thread-safe!!!
-  const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
   //
-  Int_t ncand = 0;
-  int npix = maxID-minID;
-  if (npix<1 || (ncand=gRandom->Poisson(npix*probNoisy))<1) return ncand; // decide how many noisy pixels will be added
   ncand = GenOrderedSample(npix,ncand,ordSample,ordSampleInd); 
   int* ordV = ordSample.GetArray();
   int* ordI = ordSampleInd.GetArray();
   for (int j=0;j<ncand;j++) {
-    fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],iz,ix);   // create noisy digit
-    dig.SetCoord1(iz);
-    dig.SetCoord2(ix);
-    dig.SetSignal(1);
-    dig.SetSignalPix((Int_t)AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,base,noise));
-    for (int k=knmaxtrk;k--;) {
-      dig.SetTrack(k,-3);
-      dig.SetHit(k,-1);
-    }
-    aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix,&dig);
-    if (GetDebug(2)) AliInfo(Form("Add noisy pixel %d(%d/%d) Noise=%d",ordV[ordI[j]],iz,ix,dig.GetSignalPix()));
+    fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],col,row,iCycle);   // create noisy digit
+    iCycle = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(row,col,tof);
+    UpdateMapNoise(col,row,AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,noiseMean,noiseSig),  iCycle);
   }
   return ncand;
 }
 
+
 //______________________________________________________________________
-void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old, Int_t ntrack, Int_t idhit
+void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old) 
 {
   //  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.
+  //  Note pList is goten via GetMap() and chip 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:
-  // old                  existing AliITSUSDigit
-  // Int_t ntrack         track incex number
-  // Int_t idhit          hit index number
   UInt_t col,row;
+  Int_t iCycle;
   Double_t pulse1,pulse2;
   Double_t couplR=0.0,couplC=0.0;
   Double_t xr=0.;
   //
-  fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
+  fSensMap->GetMapIndex(old->GetUniqueID(),col,row,iCycle);
+  int cycle = iCycle;
   fSimuParam->GetPixCouplingParam(couplC,couplR);
-  if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
-                               col,row,ntrack,idhit,couplC,couplR));
+  if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,couplC=%e couplR=%e",
+                               col,row,couplC,couplR));
   pulse2 = pulse1 = old->GetSignal();
   if (pulse1<fSimuParam->GetPixMinElToAdd()) return; // too small signal
   for (Int_t isign=-1;isign<=1;isign+=2) {
@@ -717,48 +653,38 @@ void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old, Int_t ntrack, Int_t i
     // loop in col direction
     int j1 = int(col) + isign;
     xr = gRandom->Rndm();
-    if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
+    if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,old->GetTrack(0),old->GetHit(0),pulse1,cycle);
     //
     // loop in row direction
     int j2 = int(row) + isign;
     xr = gRandom->Rndm();
-    if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
+    if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),old->GetTrack(0),old->GetHit(0),pulse2,cycle);
   } 
   //
 }
 
 //______________________________________________________________________
-void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, Int_t ntrack,Int_t idhit
+void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old) 
 {
   //  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:
   // old            existing AliITSUSDigit
   // ntrack         track incex number
   // idhit          hit index number
-  // module         module number
+  // chip         chip number
   //
   UInt_t col,row;
+  Int_t cycle;
+  Int_t modId = fChip->GetIndex();
   Double_t pulse1,pulse2;
   Double_t couplR=0.0,couplC=0.0;
   //
-  fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
+  fSensMap->GetMapIndex(old->GetUniqueID(),col,row,cycle);
   fSimuParam->GetPixCouplingParam(couplC,couplR);
-  if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
-                               col,row,ntrack,idhit,couplC,couplR));
+  if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,roCycle=%d)  couplC=%e couplR=%e",col,row,cycle,couplC,couplR));
  //
  if (old->GetSignal()<fSimuParam->GetPixMinElToAdd()) return; // too small signal
  for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
@@ -766,26 +692,132 @@ void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, Int_t ntrack,Int_t
    //
    int j1 = int(col)+isign;
    pulse1 *= couplC;    
-   if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(fModule))) pulse1 = old->GetSignal();
-   else UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
+   if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(modId))) pulse1 = old->GetSignal();
+   else UpdateMapSignal(UInt_t(j1),row,old->GetTrack(0),old->GetHit(0),pulse1,cycle);
    
    // loop in row direction
    int j2 = int(row) + isign;
    pulse2 *= couplR;
-   if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(fModule))) pulse2 = old->GetSignal();
-   else UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
+   if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(modId))) pulse2 = old->GetSignal();
+   else UpdateMapSignal(col,UInt_t(j2),old->GetTrack(0),old->GetHit(0),pulse2,cycle);
  } // for isign
 }
 
 //______________________________________________________________________
-void AliITSUSimulationPix::GenerateStrobePhase()
+void AliITSUSimulationPix::SetResponseParam(AliITSUParamList* resp)
+{
+  // attach response parameterisation data
+  fResponseParam = resp;
+  //
+  int spreadID = Nint(fResponseParam->GetParameter(AliITSUSimulationPix::kChargeSpreadType));
+  const char* hname = 0;
+  fSpread2DHisto = 0;
+  //
+  switch (spreadID) {
+    //
+  case kSpreadFunHisto: 
+    fSpreadFun = &AliITSUSimulationPix::SpreadFrom2DHisto;
+    hname = fResponseParam->GetParName(AliITSUSimulationPix::kChargeSpreadType);
+    if (!(fSpread2DHisto=(TH2*)fResponseParam->GetParamObject(hname))) 
+      AliFatal(Form("Did not find 2D histo %s for charge spread parameterization",hname));
+    break;
+    //
+  case kSpreadFunDoubleGauss2D: 
+    fSpreadFun = &AliITSUSimulationPix::SpreadFunDoubleGauss2D; 
+    break;
+    //
+  case kSpreadFunGauss2D: 
+    fSpreadFun = &AliITSUSimulationPix::SpreadFunGauss2D;       
+    break;
+    //
+  default: AliFatal(Form("Did not find requested spread function id=%d",spreadID));
+  }
+  //
+  int readoutType = Nint(fResponseParam->GetParameter(kReadOutSchemeType));
+  switch (readoutType) {
+  case kReadOutStrobe: 
+    fROTimeFun = &AliITSUSimulationPix::GetReadOutCycle;
+    break;
+  case kReadOutRollingShutter: 
+    fROTimeFun = &AliITSUSimulationPix::GetReadOutCycleRollingShutter;
+    break;
+  default: AliFatal(Form("Did not find requested readout time type id=%d",readoutType));
+  }
+  //___ Set the Rolling Shutter read-out window 
+  fReadOutCycleLength = fResponseParam->GetParameter(kReadOutCycleLength);
+  //___ Pixel discrimination threshold, and the S/N cut
+  fSimuParam->SetPixThreshold(fResponseParam->GetParameter(kPixNoiseMPV) *fResponseParam->GetParameter(kPixSNDisrcCut) , fResponseParam->GetParameter(kPixSNDisrcCut),-1); //for all chips
+  //___ Minimum number of electrons to add 
+  fSimuParam->SetPixMinElToAdd(fResponseParam->GetParameter(kPixMinElToAdd));
+  //___ Set the Pixel Noise MPV and Sigma (the noise distribution is Landau not Gauss due to RTN)
+  fSimuParam->SetPixNoise( fResponseParam->GetParameter(kPixNoiseMPV), fResponseParam->GetParameter(kPixNoiseSigma), -1); //for all chips
+  //___ Pixel fake hit rate
+  fSimuParam->SetPixFakeRate( fResponseParam->GetParameter(kPixFakeRate) );
+  //___ To apply the noise or not
+  if (  fResponseParam->GetParameter(kPixNoiseIsOn) > 0.01)  fSimuParam->SetPixAddNoisyFlag(kTRUE);
+  else fSimuParam->SetPixAddNoisyFlag(kFALSE);
+  //
+  if(fResponseParam->GetParameter(kPixNoiseInAllMod) > 0.01 ) fSimuParam->SetPixNoiseInAllMod(kTRUE);
+  else fSimuParam->SetPixNoiseInAllMod(kFALSE);
+  //
+  //  Double_t vGeVToQ = fSimuParam->GetGeVToCharge();
+  fGlobalChargeScale = fResponseParam->GetParameter(kSpreadFunGlobalQScale);
+    
+  AliDebug(10,Form("=============== Setting the response start ============================"));
+  AliDebug(10,Form("=============== RO type: %d",readoutType));
+  AliDebug(10,Form("=============== RO cycle lenght: %lf",fReadOutCycleLength));
+  AliDebug(10,Form("=============== Noise MPV: %lf",fResponseParam->GetParameter(kPixNoiseMPV)));
+  AliDebug(10,Form("=============== Noise Sigma: %lf",fResponseParam->GetParameter(kPixNoiseSigma)));
+  AliDebug(10,Form("=============== Fake rate: %lf",fResponseParam->GetParameter(kPixFakeRate)));
+  AliDebug(10,Form("=============== Noise On/Off: %d",fSimuParam->GetPixAddNoisyFlag()));
+  AliDebug(10,Form("=============== Noise in all mod on/off: %d",fSimuParam->GetPixNoiseInAllMod()));
+  AliDebug(10,Form("=============== Global Charge scale: %lf",fGlobalChargeScale));
+  AliDebug(10,Form("=============== Setting the response done  ============================"));
+  
+  
+  
+}
+
+//______________________________________________________________________
+Int_t AliITSUSimulationPix::GetReadOutCycleRollingShutter(Int_t row, Int_t col, Double_t hitTime)
 {
-  // Generate randomly the strobe
-  // phase w.r.t to the LHC clock
-  // Done once per event
-  const Double_t kBunchLenght = 25e-9; // LHC clock
-  fStrobePhase = ((Double_t)gRandom->Integer(fStrobeLenght))*kBunchLenght-
-    (Double_t)fStrobeLenght*kBunchLenght+
-    kBunchLenght/2;
+  //
+  // Get the read-out cycle of the hit in the given column/row of the sensor.
+  // hitTime is the time of the subhit (hit is divided to nstep charge deposit) in seconds
+  // globalPhaseShift gives the start of the RO for the cycle in pixel wrt the LHC clock
+  // GetRollingShutterWindow give the with of the rolling shutter read out window
+  //
+  double tmin = fReadOutCycleOffset + fReadOutCycleLength*(double(row)/fSeg->Npx()-1.);
+  int cycle = Nint( (hitTime-tmin)/fReadOutCycleLength - 0.5 );
+  AliDebug(3,Form("Rolling shutter at row%d/col%d: particle time: %e, tmin: %e : tmax %e -> cycle:%d",row,col,hitTime,tmin,
+                 tmin+fReadOutCycleLength,cycle));
+  return cycle;
+  //  
 }
 
+//______________________________________________________________________
+Int_t AliITSUSimulationPix::GetReadOutCycle(Int_t row, Int_t col, Double_t hitTime)
+{
+  //
+  // Check whether the hit is in the read out window of the given column/row of the sensor
+  //
+  AliDebug(3,Form("Strobe readout: row%d/col%d: particle time: %e, tmin: %e, tmax %e",
+                 row,col,hitTime,fReadOutCycleOffset,fReadOutCycleOffset+fReadOutCycleLength));
+  hitTime -= fReadOutCycleOffset+0.5*fReadOutCycleLength;
+  return (hitTime<0 || hitTime>fReadOutCycleLength) ? kMaxROCycleAccept+1 : 0;
+  //  
+}
+
+//_______________________________________________________________________
+void AliITSUSimulationPix::CalcDiodeShiftInPixel(Int_t xrow, Int_t zcol, Float_t &x, Float_t &z)
+{
+  //
+  // Calculates the shift of the diode wrt the geometric center of the pixel.
+  // It is needed for staggerred pixel layout or double diode pixels with assymetric center
+  // The shift can depend on the column or line or both...
+  // The x and z are passed in cm 
+  //
+  ((AliITSUSegmentationPix*)fSeg)->GetDiodShift(xrow,zcol,x,z);  
+  //
+}
+//_______________________________________________________________________