]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ITS/AliITSsimulationSDD.cxx
fix in calling of gaussian spread function
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSDD.cxx
index 7dd42fb91e3086ca27f02eb321761d417833046a..165efd78deef9a9bf3a24cd1c973494abf68d65f 100644 (file)
@@ -18,7 +18,7 @@
 #include <Riostream.h>
 #include <stdlib.h>
 #include <stdio.h>
-#include <string.h>
+#include <cstring>
 
 #include <TCanvas.h>
 #include <TF1.h>
 #include <TROOT.h>
 #include "AliITS.h"
 #include "AliITSMapA2.h"
-#include "AliITSRawData.h"
 #include "AliITSdigitSPD.h"
 #include "AliITSetfSDD.h"
 #include "AliITSmodule.h"
+#include "AliITShit.h"
 #include "AliITSpList.h"
 #include "AliITSCalibrationSDD.h"
-#include "AliITSsegmentationSDD.h"
+#include "AliITSresponseSDD.h"
 #include "AliITSsimulationSDD.h"
 #include "AliLog.h"
 #include "AliRun.h"
 
+using std::endl;
+using std::cout;
 ClassImp(AliITSsimulationSDD)
 ////////////////////////////////////////////////////////////////////////
 // Version: 0                                                         //
@@ -69,54 +71,11 @@ fNofMaps(0),
 fMaxNofSamples(0),
 fScaleSize(0){
     // Default constructor
-    SetScaleFourier();
     SetPerpendTracksFlag();
     SetCrosstalkFlag();
     SetDoFFT();
 }
-//______________________________________________________________________
-AliITSsimulationSDD::AliITSsimulationSDD(const AliITSsimulationSDD &source) :
-    AliITSsimulation(source),
-fITS(source.fITS),
-fHitMap2(source.fHitMap2),
-fHitSigMap2(source.fHitSigMap2),
-fHitNoiMap2(source.fHitNoiMap2),
-fElectronics(source.fElectronics),
-fInZR(source.fInZR),
-fInZI(source.fInZI),
-fOutZR(source.fOutZR),
-fOutZI(source.fOutZI),
-fAnodeFire(source.fAnodeFire),
-fHis(source.fHis),
-fFlag(source.fFlag),
-fCrosstalkFlag(source.fCrosstalkFlag),
-fDoFFT(source.fDoFFT),
-fNofMaps(source.fNofMaps),
-fMaxNofSamples(source.fMaxNofSamples),
-fScaleSize(source.fScaleSize){
-    // Copy constructor to satify Coding roules only.
-
-}
-//______________________________________________________________________
-AliITSsimulationSDD& AliITSsimulationSDD::operator=(const AliITSsimulationSDD &src){
-    // Assignment operator to satify Coding roules only.
-
-    if(this==&src) return *this;
-    Error("AliITSsimulationSDD","Not allowed to make a = with "
-          "AliITSsimulationSDD Using default creater instead");
-    return *this ;
-}
-/*
-//______________________________________________________________________
-AliITSsimulation& AliITSsimulationSDD::operator=(const AliITSsimulation &src){
-    // Assignment operator to satify Coding roules only.
 
-    if(this==&src) return *this;
-    Error("AliITSsimulationSSD","Not allowed to make a = with "
-          "AliITSsimulationSDD Using default creater instead");
-    return *this ;
-}
-*/
 //______________________________________________________________________
 AliITSsimulationSDD::AliITSsimulationSDD(AliITSDetTypeSim* dettyp):
 AliITSsimulation(dettyp),
@@ -144,13 +103,12 @@ fScaleSize(0){
 void AliITSsimulationSDD::Init(){
     // Standard Constructor
 
-    SetScaleFourier();
+    AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
+    fScaleSize = ScaleFourier(seg);
     SetPerpendTracksFlag();
     SetCrosstalkFlag();
     SetDoFFT();
 
-    AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
-    
     AliITSSimuParam* simpar = fDetType->GetSimuParam();
     fpList = new AliITSpList( seg->Npz(),
                               fScaleSize*seg->Npx() );
@@ -167,9 +125,8 @@ void AliITSsimulationSDD::Init(){
     Double_t timeStep  = (Double_t)seg->Dpx(0);
 
     if(anodePitch*(fNofMaps/2) > sddWidth) {
-        Warning("AliITSsimulationSDD",
-                "Too many anodes %d or too big pitch %f \n",
-                fNofMaps/2,anodePitch);
+      AliWarning(Form("Too many anodes %d or too big pitch %f ",
+                fNofMaps/2,anodePitch));
     } // end if
 
 
@@ -251,7 +208,6 @@ void AliITSsimulationSDD::FastFourierTransform(Double_t *real,
         m /= 2;
         m2 += m2;
     } // end for i
-  
     for(j=0; j<samples; j++) {
         Int_t j1 = j;
         Int_t p = 0;
@@ -385,246 +341,245 @@ void AliITSsimulationSDD::HitsToAnalogDigits( AliITSmodule *mod ) {
   AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
   AliITSSimuParam* simpar = fDetType->GetSimuParam();
-
   TObjArray *hits     = mod->GetHits();
-    Int_t      nhits    = hits->GetEntriesFast();
-
-    //    Int_t      arg[6]   = {0,0,0,0,0,0};
-    Int_t     nofAnodes  = fNofMaps/2;
-    Double_t  sddLength  = seg->Dx();
-    Double_t  sddWidth   = seg->Dz();
-    Double_t  anodePitch = seg->Dpz(0);
-    Double_t  timeStep   = seg->Dpx(0);
-    Double_t  driftSpeed ;  // drift velocity (anode dependent)
-    Double_t  norm       = simpar->GetSDDMaxAdc()/simpar->GetSDDDynamicRange(); //   maxadc/topValue;
-    Double_t  cHloss     = simpar->GetSDDChargeLoss();
-    Float_t   dfCoeff, s1; 
-    simpar->GetSDDDiffCoeff(dfCoeff,s1); // Signal 2d Shape
-    Double_t  eVpairs    = simpar->GetGeVToCharge()*1.0E9; // 3.6 eV by def.
-    Double_t  nsigma     = simpar->GetNSigmaIntegration(); //
-    Int_t     nlookups   = simpar->GetGausNLookUp();       //
-    Float_t   jitter     = simpar->GetSDDJitterError(); // 
+  Int_t      nhits    = hits->GetEntriesFast();
+
+  //    Int_t      arg[6]   = {0,0,0,0,0,0};
+  Int_t     nofAnodes  = fNofMaps/2;
+  Double_t  sddLength  = seg->Dx();
+  Double_t  anodePitch = seg->Dpz(0);
+  Double_t  timeStep   = seg->Dpx(0);
+  Double_t  driftSpeed ;  // drift velocity (anode dependent)
+  Double_t  nanoampToADC       = simpar->GetSDDMaxAdc()/simpar->GetSDDDynamicRange(); //   maxadc/topValue;
+  Double_t  cHloss     = simpar->GetSDDChargeLoss();
+  Float_t   dfCoeff, s1; 
+  simpar->GetSDDDiffCoeff(dfCoeff,s1); // Signal 2d Shape
+  Double_t  eVpairs    = simpar->GetGeVToCharge()*1.0E9; // 3.6 eV by def.
+  Double_t  nsigma     = simpar->GetNSigmaIntegration(); //
+  Int_t     nlookups   = simpar->GetGausNLookUp();       //
+  Float_t   jitter     = simpar->GetSDDJitterError(); // 
+  Float_t   mapsmear   = simpar->GetSDDCorrMapPrecision(); // 
+  Float_t   trigDelay  = simpar->GetSDDTrigDelay(); // compensation for MC time zero
+  if(res->IsAMAt20MHz()) trigDelay+=12.5; // compensation for discretization step
+
+  Float_t   timeZero=fDetType->GetResponseSDD()->GetTimeZero(fModule);
+  Float_t   adcscale   = fDetType->GetResponseSDD()->GetADCtokeV(fModule);
+  adcscale/=simpar->GetSDDkeVtoADC();
+
+  // Piergiorgio's part (apart for few variables which I made float
+  // when i thought that can be done
+  // Fill detector maps with GEANT hits
+  // loop over hits in the module
+  
+  const Float_t kconv = 1.0e+6;  // GeV->KeV
+  Int_t     itrack      = 0;
+  Int_t     iWing;       // which detector wing/side.
+  Int_t     ii,kk,ka,kt; // loop indexs
+  Int_t     ia,it,index; // sub-pixel integration indexies
+  Int_t     iAnode;      // anode number.
+  Int_t     timeSample;  // time buckett.
+  Int_t     anodeWindow; // anode direction charge integration width
+  Int_t     timeWindow;  // time direction charge integration width
+  Int_t     jamin,jamax; // anode charge integration window
+  Int_t     jtmin,jtmax; // time charge integration window
+  Int_t     nsplitAn;    // the number of splits in anode and time windows
+  Int_t     nsplitTb;    // the number of splits in anode and time windows
+  Int_t     nOfSplits;   // number of times track length is split into
+  Float_t   nOfSplitsF;  // Floating point version of nOfSplits.
+  Float_t   kkF;         // Floating point version of loop index kk.
+  Double_t  pathInSDD; // Track length in SDD.
+  Double_t  drPath; // average position of track in detector. in microns
+  Double_t  drTime; // Drift time
+  Double_t  avDrft;  // x position of path length segment in cm.
+  Double_t  avAnode; // Anode for path length segment in Anode number (float)
+  Double_t  zAnode;  // Floating point anode number.
+  Double_t  driftPath; // avDrft in microns.
+  Double_t  width;     // width of signal at anodes.
+  Double_t  depEnergy; // Energy deposited in this GEANT step.
+  Double_t  xL[3],dxL[3]; // local hit coordinates and diff.
+  Double_t  sigA; // sigma of signal at anode.
+  Double_t  sigT; // sigma in time/drift direction for track segment
+  Double_t  aStep,aConst; // sub-pixel size and offset anode
+  Double_t  tStep,tConst; // sub-pixel size and offset time
+  Double_t  amplitude; // signal amplitude for track segment in nanoAmpere
+  Double_t  chargeloss; // charge loss for track segment.
+  Double_t  anodeAmplitude; // signal amplitude in anode direction
+  Double_t  aExpo;          // exponent of Gaussian anode direction
+  Double_t  timeAmplitude;  // signal amplitude in time direction
+  Double_t  tExpo;          // exponent of Gaussian time direction
+  Double_t  tof;            // Time of flight in ns of this step.    
+  
+  for(ii=0; ii<nhits; ii++) {
+    if(!mod->LineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2],
+                         depEnergy,itrack)) continue;
+    Float_t xloc=xL[0];
+    if(xloc>0) iWing=0; // left side, carlos channel 0
+    else iWing=1; // right side
     
-    // Piergiorgio's part (apart for few variables which I made float
-    // when i thought that can be done
-    // Fill detector maps with GEANT hits
-    // loop over hits in the module
-
-    const Float_t kconv = 1.0e+6;  // GeV->KeV
-    Int_t     itrack      = 0;
-    Int_t     iWing;       // which detector wing/side.
-    Int_t     ii,kk,ka,kt; // loop indexs
-    Int_t     ia,it,index; // sub-pixel integration indexies
-    Int_t     iAnode;      // anode number.
-    Int_t     timeSample;  // time buckett.
-    Int_t     anodeWindow; // anode direction charge integration width
-    Int_t     timeWindow;  // time direction charge integration width
-    Int_t     jamin,jamax; // anode charge integration window
-    Int_t     jtmin,jtmax; // time charge integration window
-    Int_t     ndiv;        // Anode window division factor.
-    Int_t     nsplit;      // the number of splits in anode and time windows==1.
-    Int_t     nOfSplits;   // number of times track length is split into
-    Float_t   nOfSplitsF;  // Floating point version of nOfSplits.
-    Float_t   kkF;         // Floating point version of loop index kk.
-    Double_t  pathInSDD; // Track length in SDD.
-    Double_t  drPath; // average position of track in detector. in microns
-    Double_t  drTime; // Drift time
-    Double_t  nmul;   // drift time window multiplication factor.
-    Double_t  avDrft;  // x position of path length segment in cm.
-    Double_t  avAnode; // Anode for path length segment in Anode number (float)
-    Double_t  zAnode;  // Floating point anode number.
-    Double_t  driftPath; // avDrft in microns.
-    Double_t  width;     // width of signal at anodes.
-    Double_t  depEnergy; // Energy deposited in this GEANT step.
-    Double_t  xL[3],dxL[3]; // local hit coordinates and diff.
-    Double_t  sigA; // sigma of signal at anode.
-    Double_t  sigT; // sigma in time/drift direction for track segment
-    Double_t  aStep,aConst; // sub-pixel size and offset anode
-    Double_t  tStep,tConst; // sub-pixel size and offset time
-    Double_t  amplitude; // signal amplitude for track segment in nanoAmpere
-    Double_t  chargeloss; // charge loss for track segment.
-    Double_t  anodeAmplitude; // signal amplitude in anode direction
-    Double_t  aExpo;          // exponent of Gaussian anode direction
-    Double_t  timeAmplitude;  // signal amplitude in time direction
-    Double_t  tExpo;          // exponent of Gaussian time direction
-    //  Double_t tof;            // Time of flight in ns of this step.    
-
-    for(ii=0; ii<nhits; ii++) {
-      if(!mod->LineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2],
-                              depEnergy,itrack)) continue;
-      Float_t xloc=xL[0];
-      if(xloc>0) iWing=0; // left side, carlos channel 0
-      else iWing=1; // right side
+    Float_t zloc=xL[2]+0.5*dxL[2];
+    zAnode=seg->GetAnodeFromLocal(xloc,zloc); // anode number in the range 0.-511.
+    driftSpeed = res->GetDriftSpeedAtAnode(zAnode);
+    driftSpeed+= fDetType->GetResponseSDD()->GetDeltaVDrift(fModule,zAnode>255);
 
-      Float_t zloc=xL[2]+0.5*dxL[2];
-      zAnode=seg->GetAnodeFromLocal(xloc,zloc); // anode number in the range 0.-511.
-      driftSpeed = res->GetDriftSpeedAtAnode(zAnode);
-      if(timeStep*fMaxNofSamples < sddLength/driftSpeed) {
-       AliWarning("Time Interval > Allowed Time Interval\n");
-      }
-      depEnergy  *= kconv;
-        
-      // scale path to simulate a perpendicular track
+    if(timeStep*fMaxNofSamples < sddLength/driftSpeed) {
+      AliWarning("Time Interval > Allowed Time Interval");
+    }
+    depEnergy  *= kconv;
+    if (!depEnergy) {
+      AliDebug(1,
+              Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!",
+                   itrack,ii,mod->GetIndex()));
+      continue;
       // continue if the particle did not lose energy
       // passing through detector
-      if (!depEnergy) {
-       AliDebug(1,
-                Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!",
-                     itrack,ii,mod->GetIndex()));
-       continue;
-      } // end if !depEnergy
-      
-      xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); //
-      pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]);
+    } // end if !depEnergy
+     
+    tof=0.;
+    AliITShit* h=(AliITShit*)hits->At(ii);
+    if(h){ 
+      tof=h->GetTOF()*1E9; 
+      AliDebug(1,Form("TOF for hit %d on mod %d (particle %d)=%g",ii,fModule,h->Track(),tof));
+    }
 
-      if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); }
-      drPath = TMath::Abs(10000.*(dxL[0]+2.*xL[0])*0.5);
-      drPath = sddLength-drPath;
-      if(drPath < 0) {
+    Float_t corrx=0, corrz=0;
+    res->GetShiftsForSimulation(xL[2],xL[0],corrz,corrx,seg);
+    xL[2]-=corrz;
+    xL[0]-=corrx;
+    xL[0] += 0.0001*gRandom->Gaus( 0, mapsmear); //
+    xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); //
+
+    pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]);
+    
+    if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); }
+    drPath = TMath::Abs(10000.*(dxL[0]+2.*xL[0])*0.5);
+    drPath = sddLength-drPath;
+    if(drPath < 0) {
+      AliInfo( // this should be fixed at geometry level
+              Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e",
+                   drPath,sddLength,dxL[0],xL[0]));
+      continue;
+    } // end if drPath < 0
+    
+    // Compute number of segments to brake step path into
+    drTime = drPath/driftSpeed;  //   Drift Time
+    sigA   = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes
+    // calcuate the number of time the path length should be split into.
+    nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA);
+    if(fFlag) nOfSplits = 1;
+    
+    // loop over path segments, init. some variables.
+    depEnergy /= nOfSplits;
+    nOfSplitsF = (Float_t) nOfSplits;
+    Float_t theAverage=0.,theSteps=0.;
+    for(kk=0;kk<nOfSplits;kk++) { // loop over path segments
+      kkF       = (Float_t) kk + 0.5;
+      avDrft    = xL[0]+dxL[0]*kkF/nOfSplitsF;
+      avAnode   = xL[2]+dxL[2]*kkF/nOfSplitsF;
+      theSteps+=1.;
+      theAverage+=avAnode;
+      zAnode = seg->GetAnodeFromLocal(avDrft,avAnode);
+      driftSpeed = res->GetDriftSpeedAtAnode(zAnode);  
+      driftSpeed+= fDetType->GetResponseSDD()->GetDeltaVDrift(fModule,zAnode>255);
+      driftPath = TMath::Abs(10000.*avDrft);
+      driftPath = sddLength-driftPath;
+      if(driftPath < 0) {
        AliDebug(1, // this should be fixed at geometry level
-                Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e",
-                  drPath,sddLength,dxL[0],xL[0]));
+                Form("negative drift path driftPath=%e sddLength=%e avDrft=%e dxL[0]=%e xL[0]=%e",
+                     driftPath,sddLength,avDrft,dxL[0],xL[0]));
        continue;
-      } // end if drPath < 0
-
-        // Compute number of segments to brake step path into
-      drTime = drPath/driftSpeed;  //   Drift Time
-      sigA   = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes
-      // calcuate the number of time the path length should be split into.
-      nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA);
-      if(fFlag) nOfSplits = 1;
+      } // end if driftPath < 0
+      drTime     = driftPath/driftSpeed; // drift time for segment.
+      // Sigma along the anodes for track segment.
+      sigA       = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
+      sigT       = sigA/driftSpeed;
+
+      drTime+=tof; // take into account Time Of Flight from production point
+      drTime-=trigDelay;
+      drTime+=timeZero;
+      timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1.001); // time bin in range 1-256 !!!
+      if(zAnode>nofAnodes) zAnode-=nofAnodes;  // to have the anode number between 0. and 256.
+      iAnode = (Int_t) (1.001+zAnode); // iAnode in range 1-256 !!!!
       
-      // loop over path segments, init. some variables.
-      depEnergy /= nOfSplits;
-      nOfSplitsF = (Float_t) nOfSplits;
-      Float_t theAverage=0.,theSteps=0.;
-      for(kk=0;kk<nOfSplits;kk++) { // loop over path segments
-       kkF       = (Float_t) kk + 0.5;
-       avDrft    = xL[0]+dxL[0]*kkF/nOfSplitsF;
-       avAnode   = xL[2]+dxL[2]*kkF/nOfSplitsF;
-       theSteps+=1.;
-       theAverage+=avAnode;
-       zAnode = seg->GetAnodeFromLocal(avDrft,avAnode);
-       driftSpeed = res->GetDriftSpeedAtAnode(zAnode); 
-       driftPath = TMath::Abs(10000.*avDrft);
-       driftPath = sddLength-driftPath;
-       if(driftPath < 0) {
-         AliDebug(1, // this should be fixed at geometry level
-                  Form("negative drift path driftPath=%e sddLength=%e avDrft=%e dxL[0]=%e xL[0]=%e",
-                       driftPath,sddLength,avDrft,dxL[0],xL[0]));
-         continue;
-       } // end if driftPath < 0
-       drTime     = driftPath/driftSpeed; // drift time for segment.
-       timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1); // time bin in range 1-256 !!!
-       if(timeSample > fScaleSize*fMaxNofSamples) {
-         AliWarning(Form("Wrong Time Sample: %e",timeSample));
-         continue;
-       } // end if timeSample > fScaleSize*fMaxNoofSamples
-       
-       if(zAnode>nofAnodes) zAnode-=nofAnodes;  // to have the anode number between 0. and 256.
-       if(zAnode*anodePitch > sddWidth || zAnode*anodePitch < 0.) 
-         AliWarning(Form("Exceeding sddWidth=%e Z = %e",sddWidth,zAnode*anodePitch));
-       iAnode = (Int_t) (1.+zAnode); // iAnode in range 1-256 !!!!
-       if(iAnode < 1 || iAnode > nofAnodes) {
-         AliWarning(Form("Wrong iAnode: 1<%d>%d  (xanode=%e)",iAnode,nofAnodes, zAnode));
-         continue;
-       } // end if iAnode < 1 || iAnode > nofAnodes
-
-       // store straight away the particle position in the array
-       // of particles and take idhit=ii only when part is entering (this
-       // requires FillModules() in the macro for analysis) :
-       
-       // Sigma along the anodes for track segment.
-       sigA       = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
-       sigT       = sigA/driftSpeed;
        // Peak amplitude in nanoAmpere
-       amplitude  = fScaleSize*160.*depEnergy/
-         (timeStep*eVpairs*2.*acos(-1.)*sigT*sigA);
-       amplitude *= timeStep/25.; // WARNING!!!!! Amplitude scaling to 
-       // account for clock variations 
-       // (reference value: 40 MHz)
-       chargeloss = 1.-cHloss*driftPath/1000.;
-       amplitude *= chargeloss;
-       width  = 2.*nsigma/(nlookups-1);
-       // Spread the charge 
-       // Pixel index
-       ndiv = 2;
-       nmul = 3.; 
-       if(drTime > 1200.) { 
-         ndiv = 4;
-         nmul = 1.5;
-       } // end if drTime > 1200.
-       // Sub-pixel index
-       nsplit = 4; // hard-wired //nsplit=4;nsplit = (nsplit+1)/2*2;
-       // Sub-pixel size see computation of aExpo and tExpo.
-       aStep  = anodePitch/(nsplit*fScaleSize*sigA);
-       aConst = zAnode*anodePitch/sigA;
-       tStep  = timeStep/(nsplit*fScaleSize*sigT);
-       tConst = drTime/sigT;
-       // Define SDD window corresponding to the hit
-       anodeWindow = (Int_t)(fScaleSize*nsigma*sigA/anodePitch+1);
-       timeWindow  = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
-       jamin = (iAnode - anodeWindow/ndiv - 2)*fScaleSize*nsplit +1;
-       jamax = (iAnode + anodeWindow/ndiv + 1)*fScaleSize*nsplit;
-       if(jamin <= 0) jamin = 1;
-       if(jamax > fScaleSize*nofAnodes*nsplit) 
-         jamax = fScaleSize*nofAnodes*nsplit;
-       // jtmin and jtmax are Hard-wired
-       jtmin = (Int_t)(timeSample-timeWindow*nmul-1)*nsplit+1;
-       jtmax = (Int_t)(timeSample+timeWindow*nmul)*nsplit;
-       if(jtmin <= 0) jtmin = 1;
-       if(jtmax > fScaleSize*fMaxNofSamples*nsplit) 
-         jtmax = fScaleSize*fMaxNofSamples*nsplit;
-       // Spread the charge in the anode-time window
-       for(ka=jamin; ka <=jamax; ka++) {
-         ia = (ka-1)/(fScaleSize*nsplit) + 1;
-         if(ia <= 0) {
-           Warning("HitsToAnalogDigits","ia < 1: ");
-           continue;
-         } // end if
-         if(ia > nofAnodes) ia = nofAnodes;
-         aExpo     = (aStep*(ka-0.5)-aConst);
-         if(TMath::Abs(aExpo) > nsigma)  anodeAmplitude = 0.;
-         else {
-           Int_t theBin = (Int_t) ((aExpo+nsigma)/width+0.5);
-           anodeAmplitude = amplitude*simpar->GetGausLookUp(theBin);
-         } // end if TMath::Abs(aEspo) > nsigma
-         // index starts from 0
-         index = iWing*nofAnodes+ia-1;
-         if(anodeAmplitude){
-           for(kt=jtmin; kt<=jtmax; kt++) {
-             it = (kt-1)/nsplit+1;  // it starts from 1
-             if(it<=0){
-               Warning("HitsToAnalogDigits","it < 1:");
-               continue;
-             } // end if 
-             if(it>fScaleSize*fMaxNofSamples)
-               it = fScaleSize*fMaxNofSamples;
-             tExpo    = (tStep*(kt-0.5)-tConst);
-             if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.;
-             else {
-               Int_t theBin = (Int_t) ((tExpo+nsigma)/width+0.5);
-               timeAmplitude = anodeAmplitude*simpar->GetGausLookUp(theBin);
-             } // end if TMath::Abs(tExpo) > nsigma
-             // build the list of Sdigits for this module        
-             //                    arg[0]     = index;
-             //                    arg[1]     = it;
-             //                    arg[2]     = itrack; // track number
-             //                    arg[3]     = ii-1; // hit number.
-             timeAmplitude *= norm;
-             timeAmplitude *= 10;
-             //         ListOfFiredCells(arg,timeAmplitude,alst,padr);
-             Double_t charge = timeAmplitude;
-             charge += fHitMap2->GetSignal(index,it-1);
-             fHitMap2->SetHit(index, it-1, charge);
-             fpList->AddSignal(index,it-1,itrack,ii-1,
-                               mod->GetIndex(),timeAmplitude);
-             fAnodeFire[index] = kTRUE;
-           }  // end loop over time in window               
-         } // end if anodeAmplitude 
-       } // loop over anodes in window
-      } // end loop over "sub-hits"
-    } // end loop over hits
+      amplitude  = fScaleSize*160.*depEnergy/
+       (timeStep*eVpairs*2.*acos(-1.));
+      chargeloss = 1.-cHloss*driftPath/1000.;
+      amplitude *= chargeloss;
+      amplitude *= adcscale;
+      width  = 2.*nsigma/(nlookups-1);
+      // Spread the charge 
+      nsplitAn = 4; 
+      nsplitTb=4;
+      aStep  = anodePitch/(nsplitAn*sigA);
+      aConst = zAnode*anodePitch/sigA;
+      tStep  = timeStep/(nsplitTb*fScaleSize*sigT);
+      tConst = drTime/sigT;
+      // Define SDD window corresponding to the hit
+      anodeWindow = (Int_t)(nsigma*sigA/anodePitch+1);
+      timeWindow  = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
+      jamin = (iAnode - anodeWindow - 2)*nsplitAn+1;
+      if(jamin <= 0) jamin = 1;
+      if(jamin > nofAnodes*nsplitAn){ 
+       AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode min=%d",jamin));
+       continue;
+      }
+      jamax = (iAnode + anodeWindow + 2)*nsplitAn;
+      if(jamax > nofAnodes*nsplitAn) jamax = nofAnodes*nsplitAn;
+      if(jamax <=0){ 
+       AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode max=%d",jamax));
+       continue;
+      }
+      jtmin = (Int_t)(timeSample-timeWindow-2)*nsplitTb+1;
+      if(jtmin <= 0) jtmin = 1;
+      if(jtmin > fScaleSize*fMaxNofSamples*nsplitTb){ 
+       AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample min=%d  tof=%f",jtmin,tof));
+       continue; 
+      }
+      jtmax = (Int_t)(timeSample+timeWindow+2)*nsplitTb;
+      if(jtmax > fScaleSize*fMaxNofSamples*nsplitTb) jtmax = fScaleSize*fMaxNofSamples*nsplitTb;
+      if(jtmax <= 0){
+       AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample max=%d  tof=%f",jtmax,tof));
+       continue; 
+      }
+
+      // Spread the charge in the anode-time window
+      for(ka=jamin; ka <=jamax; ka++) {          
+       ia = (ka-1)/nsplitAn + 1;
+       if(ia <= 0) ia=1; 
+       if(ia > nofAnodes) ia = nofAnodes;
+       aExpo     = (aStep*(ka-0.5)-aConst);
+       if(TMath::Abs(aExpo) > nsigma)  anodeAmplitude = 0.;
+       else {
+         Int_t theBin = (Int_t) ((aExpo+nsigma)/width+0.5);
+         anodeAmplitude = amplitude*simpar->GetGausLookUp(theBin);
+       }
+       // index starts from 0
+       index = iWing*nofAnodes+ia-1;
+       if(anodeAmplitude){
+         for(kt=jtmin; kt<=jtmax; kt++) {
+           it = (kt-1)/nsplitTb+1;  // it starts from 1
+           if(it<=0) it=1;
+           if(it>fScaleSize*fMaxNofSamples)
+             it = fScaleSize*fMaxNofSamples;
+           tExpo    = (tStep*(kt-0.5)-tConst);
+           if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.;
+           else {
+             Int_t theBin = (Int_t) ((tExpo+nsigma)/width+0.5);
+             timeAmplitude = anodeAmplitude*simpar->GetGausLookUp(theBin)*aStep*tStep;
+           }
+           timeAmplitude *= nanoampToADC;
+           //         ListOfFiredCells(arg,timeAmplitude,alst,padr);
+           Double_t charge = timeAmplitude;
+           charge += fHitMap2->GetSignal(index,it-1);
+           fHitMap2->SetHit(index, it-1, charge);
+           fpList->AddSignal(index,it-1,itrack,ii-1,
+                             mod->GetIndex(),timeAmplitude);
+           fAnodeFire[index] = kTRUE;
+         }  // end loop over time in window               
+       } // end if anodeAmplitude 
+      } // loop over anodes in window
+    } // end loop over "sub-hits"
+  } // end loop over hits
 }
 
 //____________________________________________
@@ -684,12 +639,16 @@ void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAdd
   Int_t i,k,kk;
   AliITSSimuParam* simpar = fDetType->GetSimuParam();
   Float_t maxadc = simpar->GetSDDMaxAdc();    
+  Int_t nGroup=fScaleSize;
+  if(res->IsAMAt20MHz()){
+    nGroup=fScaleSize/2;
+  }
 
   for (i=0;i<fNofMaps;i++) {
     if( !fAnodeFire[i] ) continue;
     baseline = res->GetBaseline(i);
     noise = res->GetNoise(i);
-    gain = res->GetChannelGain(i);
+    gain = res->GetChannelGain(i)/fDetType->GetAverageGainSDD();
     if(res->IsBad()) gain=0.;
     if( res->IsChipBad(res->GetChip(i)) )gain=0.;
     for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
@@ -701,7 +660,7 @@ void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAdd
       }
       fInZI[k]  = 0.;
     } // end for k
-    if(!fDoFFT) {
+    if(!fDoFFT) {      
       for(k=0; k<fMaxNofSamples; k++) {
        Double_t newcont = 0.;
        Double_t maxcont = 0.;
@@ -712,7 +671,7 @@ void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAdd
        newcont = maxcont;
        if (newcont >= maxadc) newcont = maxadc -1;
        if(newcont >= baseline){
-         Warning("","newcont=%d>=baseline=%d",newcont,baseline);
+         Warning("","newcont=%f>=baseline=%f",newcont,baseline);
        } // end if
          // back to analog: ?
        fHitMap2->SetHit(i,k,newcont);
@@ -729,7 +688,7 @@ void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAdd
       for(k=0; k<fMaxNofSamples; k++) {
        Double_t newcont1 = 0.;
        Double_t maxcont1 = 0.;
-       for(kk=0;kk<fScaleSize;kk++) {
+       for(kk=0;kk<nGroup;kk++) {
          newcont1 = fOutZR[fScaleSize*k+kk];
          if(newcont1 > maxcont1) maxcont1 = newcont1;
        } // end for kk
@@ -749,10 +708,6 @@ void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
   
     // create and inizialice crosstalk map
     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
-    if( ctk == NULL ) {
-        Error( "ApplyCrosstalk", "no memory for temporal map: exit \n" );
-        return;
-    }
     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
     for( Int_t z=0; z<fNofMaps; z++ ) {
@@ -781,11 +736,6 @@ void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
                         // make smooth derivative
                         Float_t* dev = new Float_t[fMaxNofSamples+1];
                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
-                        if( ctk == NULL ) {
-                            Error( "ApplyCrosstalk", 
-                                   "no memory for temporal array: exit \n" );
-                            return;
-                        }
                         for( Int_t i=tstart; i<tstop; i++ ) {   
                             if( i > 2 && i < fMaxNofSamples-2 )
                                 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 ) 
@@ -893,12 +843,12 @@ void AliITSsimulationSDD::Compress2D(){
        if(sS>tL) nLow++;
        if(sS>tH) nHigh++;
        
-       if(nLow>=3 && nHigh>=1){
+       if(nLow>=2 && nHigh>=1){
          Int_t signal=(Int_t)cC;
          Int_t signalc = Convert10to8(signal);
          Int_t signale = Convert8to10(signalc);
          signalc-=tL; // subtract low threshold after 10 to 8 bit compression
-         AddDigit(ian,itb,signalc,signale);  // store C 
+         if(signalc>=4) AddDigit(ian,itb,signalc,signale);  // store C 
        }
       }
     }
@@ -920,18 +870,16 @@ void AliITSsimulationSDD::StoreAllDigits(){
 } 
 //______________________________________________________________________
 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
-    // Creates histograms of maps for debugging
-    Int_t i;
-
-    fHis=new TObjArray(fNofMaps);
-    for (i=0;i<fNofMaps;i++) {
-        TString sddName("sdd_");
-        Char_t candNum[4];
-        sprintf(candNum,"%d",i+1);
-        sddName.Append(candNum);
-        fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
-                             0.,(Float_t) scale*fMaxNofSamples), i);
-    } // end for i
+  // Creates histograms of maps for debugging
+  Int_t i;
+  
+  fHis=new TObjArray(fNofMaps);
+  for (i=0;i<fNofMaps;i++) {
+    TString sddName;
+    sddName.Form("sdd_%d",i+1);
+    fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
+                        0.,(Float_t) scale*fMaxNofSamples), i);
+  } // end for i
 }
 //______________________________________________________________________
 void AliITSsimulationSDD::FillHistograms(){