add getter to cut on n cells
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSSD.cxx
index 27650cf..b28c288 100644 (file)
+/**************************************************************************
+ * 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 <stdio.h>
+#include <stdlib.h>
+#include <Riostream.h>
 #include <TObjArray.h>
-
-#include "AliITSsimulationSSD.h"
-#include "AliITSdictSSD.h"
-#include "AliITSdcsSSD.h"
+#include <TRandom.h>
+
+#include <TGeoGlobalMagField.h>
+#include "AliITSmodule.h"
+#include "AliITSMapA2.h"
+#include "AliITSpList.h"
+#include "AliITSCalibrationSSD.h"
+#include "AliITSsegmentationSSD.h"
+//#include "AliITSdcsSSD.h"
 #include "AliITS.h"
+#include "AliITShit.h"
+#include "AliITSdigitSSD.h"
 #include "AliRun.h"
+#include "AliMagF.h"
+#include "AliITSgeom.h"
+#include "AliITSsimulationSSD.h"
+#include "AliITSTableSSD.h"
+#include <TF1.h>
+#include "AliMathBase.h"
+
+using std::endl;
+using std::cout;
+ClassImp(AliITSsimulationSSD)
+////////////////////////////////////////////////////////////////////////
+//                                                                    //
+// Author: Enrico Fragiacomo                                          //
+//         enrico.fragiacomo@ts.infn.it                               //
+// Last revised: june 2008                                            // 
+//                                                                    //
+// AliITSsimulationSSD is the simulation of SSD.                     //
+////////////////////////////////////////////////////////////////////////
+
+//----------------------------------------------------------------------
+AliITSsimulationSSD::AliITSsimulationSSD():AliITSsimulation(),
+                                          //fDCS(0),
+fMapA2(0),
+fIonE(0.0),
+fDifConst(),
+fDriftVel(),
+fTimeResponse(NULL),
+fLorentz(kFALSE),
+fTanLorAngP(0),
+fTanLorAngN(0)
+{
+    //default Constructor
+    //Inputs:
+    // none.
+    // Outputs:
+    // none.
+    // Return:
+    //  A default construction AliITSsimulationSSD class
+}
+//----------------------------------------------------------------------
+AliITSsimulationSSD::AliITSsimulationSSD(AliITSDetTypeSim* dettyp):
+AliITSsimulation(dettyp),
+//fDCS(0),
+fMapA2(0),
+fIonE(0.0),
+fDifConst(),
+fDriftVel(),
+fTimeResponse(NULL),
+fLorentz(kFALSE),
+fTanLorAngP(0),
+fTanLorAngN(0)
+{
+    // Constructor 
+    // Input:
+    //   AliITSDetTypeSim    Pointer to the SSD dettype to be used
+    // Outputs:
+    //   none.
+    // Return
+    //   A standard constructed AliITSsimulationSSD class
+
+  fTimeResponse = new TF1("ftimeresponse",".5*x*exp(1.-.5*x)");
+    Init();
+}
+//----------------------------------------------------------------------
+void AliITSsimulationSSD::Init(){
+  // Inilizer, Inilizes all of the variable as needed in a standard place.
+  // Input:
+  //   AliITSsegmentationSSD *seg  Pointer to the SSD segmentation to be used
+  //   AliITSCalibrationSSD   *resp Pointer to the SSD responce class to be used
+  // Outputs:
+  //   none.
+  // Return
+  //   none.
+  AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2);
+  AliITSSimuParam* simpar = fDetType->GetSimuParam();
+  
+  SetDriftVelocity(); // use default values in .h file
+  SetIonizeE();       // use default values in .h file
+  SetDiffConst();     // use default values in .h file
+  fpList           = new AliITSpList(2,GetNStrips());
+  fMapA2           = new AliITSMapA2(seg);
+  SetLorentzDrift(simpar->GetSSDLorentzDrift());
+  if (fLorentz) SetTanLorAngle();
+}
 
+//______________________________________________________________________
+Bool_t AliITSsimulationSSD::SetTanLorAngle() {
+    // This function set the Tangent of the Lorentz angles. 
+    // output: Bool_t : kTRUE in case of success
+    //
 
-ClassImp(AliITSsimulationSSD);
-//------------------------------------------------------------
-AliITSsimulationSSD::AliITSsimulationSSD(AliITSsegmentation *seg,
-                                         AliITSresponse *resp){
-  // Constructor
+    if(!fDetType) {
+      AliError("AliITSsimulationSPD::SetTanLorAngle: AliITSDetTypeSim* fDetType not set ");
+      return kFALSE;}
 
-    fSegmentation = seg;
-    fResponse = resp;
-    fDCS = new AliITSdcsSSD(seg,resp); 
+    AliITSSimuParam* simpar = fDetType->GetSimuParam();
+    AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
+    if (!fld) AliFatal("The field is not initialized");
+    Double_t bz = fld->SolenoidField();
 
-    fNstrips = fSegmentation->Npx();
-    fPitch = fSegmentation->Dpx(0);
-    
-    fP = new TArrayF(fNstrips+1); 
-    fN = new TArrayF(fNstrips+1);
-     
-    fTracksP = new AliITSdictSSD[fNstrips+1];
-    fTracksN = new AliITSdictSSD[fNstrips+1];
+    fTanLorAngN = TMath::Tan( simpar->LorentzAngleElectron(bz) );
+    fTanLorAngP = TMath::Tan( simpar->LorentzAngleHole(bz) );
 
-    
-    fSteps  = 10;   // still hard-wired - set in SetDetParam and get it via  
-                     // fDCS together with the others eventually    
+    return kTRUE;
+}
 
+//______________________________________________________________________
+AliITSsimulationSSD& AliITSsimulationSSD::operator=(
+                                         const AliITSsimulationSSD &s){
+  // Operator =
+
+  if(this==&s) return *this;
+
+  //  this->fDCS         = new AliITSdcsSSD(*(s.fDCS));
+  this->fMapA2       = s.fMapA2;
+  this->fIonE        = s.fIonE;
+  this->fDifConst[0] = s.fDifConst[0];
+  this->fDifConst[1] = s.fDifConst[1];
+  this->fDriftVel[0] = s.fDriftVel[0];
+  this->fDriftVel[1] = s.fDriftVel[1];
+  this->fTimeResponse = s.fTimeResponse;
+  this->fLorentz   = s.fLorentz;
+  this->fTanLorAngP = s.fTanLorAngP;
+  this->fTanLorAngN = s.fTanLorAngN;
+  return *this;
 }
-//___________________________________________________________________________
-AliITSsimulationSSD& AliITSsimulationSSD::operator=(AliITSsimulationSSD 
-                                                                      &source){
-// Operator =
-    if(this==&source) return *this;
-
-    this->fDCS = new AliITSdcsSSD(*(source.fDCS));
-    this->fN   = new TArrayF(*(source.fN));
-    this->fP   = new TArrayF(*(source.fP));
-    this->fTracksP = new AliITSdictSSD(*(source.fTracksP));
-    this->fTracksN = new AliITSdictSSD(*(source.fTracksN));
-    this->fNstrips = source.fNstrips;
-    this->fPitch   = source.fPitch;
-    this->fSteps   = source.fSteps;
-    return *this;
+/*
+//______________________________________________________________________
+AliITSsimulation& AliITSsimulationSSD::operator=(
+                                         const AliITSsimulation &s){
+  // Operator =
+
+  if(this==&s) return *this;
+  Error("AliITSsimulationSSD","Not allowed to make a = with "
+       "AliITSsimulationSSD Using default creater instead");
+  
+  return *this;
 }
-//_____________________________________________________________
-AliITSsimulationSSD::AliITSsimulationSSD(AliITSsimulationSSD &source){
+*/
+//______________________________________________________________________
+AliITSsimulationSSD::AliITSsimulationSSD(const AliITSsimulationSSD &source):
+    AliITSsimulation(source),
+fMapA2(source.fMapA2),
+fIonE(source.fIonE),
+fDifConst(),
+fDriftVel(),
+fTimeResponse(source.fTimeResponse),
+fLorentz(source.fLorentz),
+fTanLorAngP(source.fTanLorAngP),
+fTanLorAngN(source.fTanLorAngN)
+{
   // copy constructor
-   *this = source;
+  fDifConst[0] = source.fDifConst[0];
+  fDifConst[1] = source.fDifConst[1];
+  fDriftVel[0] = source.fDriftVel[0];
+  fDriftVel[1] = source.fDriftVel[1];
 }
-//____________________________________________________________________________
+//______________________________________________________________________
 AliITSsimulationSSD::~AliITSsimulationSSD() {
-  // anihilator    
+  // destructor
+  delete fMapA2;
+  delete fTimeResponse;
+  //delete fDCS;
+}
+//______________________________________________________________________
+void AliITSsimulationSSD::InitSimulationModule(Int_t module,Int_t event){
+    // Creates maps to build the list of tracks for each sumable digit
+    // Inputs:
+    //   Int_t module    // Module number to be simulated
+    //   Int_t event     // Event number to be simulated
+    // Outputs:
+    //   none.
+    // Return
+    //    none.
+
+    SetModuleNumber(module);
+    SetEventNumber(event);
+    fMapA2->ClearMap();
+    fpList->ClearMap();
+}
+//______________________________________________________________________
+void AliITSsimulationSSD::FinishSDigitiseModule(){
+    // Does the Sdigits to Digits work
+    // Inputs:
+    //   none.
+    // Outputs:
+    //   none.
+    // Return:
+    //   none.
+
+  FillMapFrompList(fpList);  // need to check if needed here or not????
+  SDigitToDigit(fModule,fpList);
+  fpList->ClearMap();
+  fMapA2->ClearMap();
+}
+//______________________________________________________________________
+void AliITSsimulationSSD::DigitiseModule(AliITSmodule *mod,Int_t,Int_t) {
+  // Digitizes hits for one SSD module
+  SetModuleNumber(mod->GetIndex());
+  
+  HitsToAnalogDigits(mod,fpList);
+  SDigitToDigit(GetModuleNumber(),fpList);
+  
+  fpList->ClearMap();
+  fMapA2->ClearMap();
+}
+//______________________________________________________________________
+void AliITSsimulationSSD::SDigitiseModule(AliITSmodule *mod,Int_t,Int_t) {
+  // Produces Summable/Analog digits and writes them to the SDigit tree. 
 
-    if(fP) delete fP;
-    if(fN) delete fN;
-    
-    if(fTracksP) delete fTracksP;
-    if(fTracksN) delete fTracksN;
-
-    delete fDCS;
-} 
-//_______________________________________________________________
-//
-// Hit to digit
-//_______________________________________________________________
-//
-void AliITSsimulationSSD::DigitiseModule(AliITSmodule *mod,Int_t module,
-                                         Int_t dummy) {
-  // Digitizes one SSD module of hits.
+    HitsToAnalogDigits(mod,fpList);
+
+    WriteSDigits(fpList);
     
-    TObjArray *hits = mod->GetHits();
-    Int_t nhits = hits->GetEntriesFast();
-    if (!nhits) return;
-
-    Int_t i;
-    for(i=0; i<fNstrips; i++) {
-       (*fP)[i] = 0;
-       (*fN)[i] = 0;
-       fTracksP[i].ZeroTracks();
-       fTracksN[i].ZeroTracks();
+    fpList->ClearMap();
+    fMapA2->ClearMap();
+}
+//______________________________________________________________________
+void AliITSsimulationSSD::SDigitToDigit(Int_t module,AliITSpList *pList){
+  // Takes the pList and finishes the digitization.
+  
+  ApplyNoise(pList,module);
+  ApplyCoupling(pList,module);
+  ApplyDeadChannels(module);
+  
+  ChargeToSignal(module,pList);
+}
+//______________________________________________________________________
+void AliITSsimulationSSD::HitsToAnalogDigits(AliITSmodule *mod,
+                                             AliITSpList *pList){
+    // Loops over all hits to produce Analog/floating point digits. This
+    // is also the first task in producing standard digits.
+  Int_t lasttrack     = -2;
+  Int_t idtrack       = -2;
+  Double_t x0=0.0, y0=0.0, z0=0.0;
+  Double_t x1=0.0, y1=0.0, z1=0.0;
+  Double_t de=0.0;
+  Int_t module = mod->GetIndex();
+  Double_t tof = 0.;
+  
+  
+  AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2);
+  
+  TObjArray *hits = mod->GetHits();
+  Int_t nhits     = hits->GetEntriesFast();
+  if (nhits<=0) return;
+  AliITSTableSSD * tav = new AliITSTableSSD(GetNStrips());
+  module = mod->GetIndex();
+  if ( mod->GetLayer() == 6 ) seg->SetLayer(6);
+  if ( mod->GetLayer() == 5 ) seg->SetLayer(5);
+
+  for(Int_t i=0; i<nhits; i++) {    
+    // LineSegmentL returns 0 if the hit is entering
+    // If hits is exiting returns positions of entering and exiting hits
+    // Returns also energy loss
+    if(GetDebug(4)){
+      cout << i << " ";
+      cout << mod->GetHit(i)->GetXL() << " "<<mod->GetHit(i)->GetYL();
+      cout << " " << mod->GetHit(i)->GetZL();
+      cout << endl;
+    } // end if
+    if (mod->LineSegmentL(i, x0, x1, y0, y1, z0, z1, de, idtrack)) {
+
+      // Scale down dE/dx according to the hit's TOF wrt to the trigger
+      // Necessary for pileup simulation
+      // EF - 21/04/09
+      tof = mod->GetHit(i)->GetTOF();
+      tof *= 1.E+6; // convert time in microsecond
+      if(tof<2.) de = de * fTimeResponse->Eval(-1.*tof+2.);
+      else de = 0.;
+      //
+
+      HitToDigit(module, x0, y0, z0, x1, y1, z1, de,tav);
+      if (lasttrack != idtrack || i==(nhits-1)) {
+       GetList(idtrack,i,module,pList,tav);
+      } // end if
+      lasttrack=idtrack;
+    } // end if
+  }  // end loop over hits
+  delete tav; tav=0;
+  return;
+}
+//----------------------------------------------------------------------
+void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0, 
+                                     Double_t z0, Double_t x1, Double_t y1, 
+                                     Double_t z1, Double_t de,
+                                     AliITSTableSSD *tav) {
+  
+  // hit to digit conversion
+  
+  AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2);
+  // Turns hits in SSD module into one or more digits.
+  //Float_t tang[2] = {0.0,0.0};
+  //seg->Angles(tang[0], tang[1]);//stereo<<->tan(stereo)~=stereo
+  Double_t x, y, z;
+  Double_t dex=0.0, dey=0.0, dez=0.0; 
+  Double_t pairs; // pair generation energy per step.
+  Double_t sigma[2] = {0.,0.};// standard deviation of the diffusion gaussian
+  Double_t tdrift[2] = {0.,0.}; // time of drift
+  Double_t w;
+  Double_t inf[2], sup[2], par0[2];                 
+  // Set up corrections for Lorentz drift (ExB)
+  Double_t tanLorAngP = fTanLorAngP;
+  Double_t tanLorAngN = fTanLorAngN;
+  if(seg->GetLayer()==6) {
+    tanLorAngP = -1.*fTanLorAngP;
+    tanLorAngN = -1.*fTanLorAngN;
+  }
+
+  // Steps in the module are determined "manually" (i.e. No Geant)
+  // NumOfSteps divide path between entering and exiting hits in steps 
+  Int_t numOfSteps = NumOfSteps(x1, y1, z1, dex, dey, dez);
+  // Enery loss is equally distributed among steps
+  de    = de/numOfSteps;
+  pairs = de/GetIonizeE(); // e-h pairs generated
+
+  //-----------------------------------------------------
+  // stepping
+  //-----------------------------------------------------
+  for(Int_t j=0; j<numOfSteps; j++) {     // stepping
+
+    x = x0 + (j+0.5)*dex;
+    y = y0 + (j+0.5)*dey;
+    if ( y > (seg->Dy()/2+10)*1.0E-4 ) {
+      // check if particle is within the detector
+      Warning("HitToDigit",
+             "hit out of detector y0=%e,y=%e,dey=%e,j =%d module=%d,  exceed=%e",
+             y0,y,dey,j,module, y-(seg->Dy()/2+10)*1.0E-4);
+      return;
+    } // end if
+    z = z0 + (j+0.5)*dez;
+
+    if(GetDebug(4)) cout <<"HitToDigit "<<x<<" "<<y<<" "<<z<< " "
+                        <<dex<<" "<<dey<<" "<<dez<<endl;
+
+    if(seg->GetLayer()==6) {
+      y=-y; // Lay6 module has sensor up-side-down!!!
     }
     
-    for(i=0; i<nhits; i++) {
-       Int_t idtrack=mod->GetHitTrackIndex(i);  
-       HitToDigit(i,idtrack,nhits,hits);
+    Int_t k;
+    //---------------------------------------------------------
+    // Pside
+    //------------------------------------------------------------
+    k=0;
+
+    // w is the coord. perpendicular to the strips
+    //    Float_t xp=x*1.e+4,zp=z*1.e+4; // microns    
+    Float_t xp=x,zp=z; 
+
+    // correction for the Lorentz's angle
+    if(fLorentz) {
+      Float_t deltaxp = (y+(seg->Dy()*1.0E-4)/2)*tanLorAngP;
+      xp+=deltaxp;  
     }
-   
-
-    ApplyNoise();
-    ApplyCoupling();    
-    ApplyThreshold();
-    ApplyDAQ();
-      
 
-}
-
-//---------------------------------------------------------------
-
-void AliITSsimulationSSD::HitToDigit(Int_t & hitNo,Int_t idtrack,
-                                     Int_t nhits,TObjArray *hits) {
-  // Turns one or more hits in an SSD module into one or more digits.
-     
-    Int_t      stripP, stripN, i;
-    Float_t    dsP, dsN;
-    Float_t    sP, sN;
-    Float_t    eP, eN;
-    Float_t    arrayEP[10];         // hard-wired number of steps
-    Float_t    arrayEN[10];
-    Int_t      track = -1;
-       
-    Float_t    ionization = 0;
-    Float_t    signal;
+    seg->GetPadTxz(xp,zp);
     
-    AliITSdictSSD *dict;  
+    // calculate drift time
+    // y is the minimum path
+    tdrift[0] = (y+(seg->Dy()*1.0E-4)/2)/GetDriftVelocity(0);
     
-    // check if this is the right order !!!!!
-
-    AliITShit *hitI = (AliITShit*)hits->At(hitNo++);
-    AliITShit *hitE = (AliITShit*)hits->At(hitNo);
-
-
-    while (!((hitE->StatusExiting()) || 
-             (hitE->StatusDisappeared()) ||
-             (hitE->StatusStop()))) {
-       
-        if (++hitNo<nhits) {
-           ionization = hitE->GetIonization(); 
-           hitE = (AliITShit*)hits->At(hitNo);
-        }
-    }   
+    w = xp; // P side strip number
     
-        
-    if (hitI->GetTrack() == hitE->GetTrack()) 
-       //track = idtrack;
-       track = hitI->GetTrack();
-    else 
-       printf("!!! Emergency !!!\n");
-     
-            
-    ionization += hitE->GetIonization();
-       
-    const Float_t kconvm=10000.;  // cm -> microns
-
-    Float_t xI, yI, zI;
-    hitI->GetPositionL(xI, yI, zI);
+    if((w<(-0.5)) || (w>(GetNStrips()-0.5))) {
+      // this check rejects hits in regions not covered by strips
+      // 0.5 takes into account boundaries 
+      if(GetDebug(4)) cout << "Dead SSD region, x,z="<<x<<","<<z<<endl;
+      return; // There are dead region on the SSD sensitive volume!!!
+    } // end if
+    // sigma is the standard deviation of the diffusion gaussian
+    if(tdrift[k]<0) return;
     
-    xI *= kconvm;
-    yI *= kconvm;
-    zI *= kconvm;
+    sigma[k] = TMath::Sqrt(2*GetDiffConst(k)*tdrift[k]);
+    sigma[k] /= (GetStripPitch()*1.0E-4);  //units of Pitch
     
-    Float_t xE, yE, zE;
-    hitE->GetPositionL(xE, yE, zE);
+    if(sigma[k]==0.0) {        
+      Error("HitToDigit"," sigma[%d]=0",k);
+      exit(0);
+    } // end if
     
-    xE *= kconvm;
-    yE *= kconvm;
-    zE *= kconvm;
+    par0[k] = pairs;
+    // we integrate the diffusion gaussian from -3sigma to 3sigma 
+    inf[k] = w - 3*sigma[k]; // 3 sigma from the gaussian average  
+    sup[k] = w + 3*sigma[k]; // 3 sigma from the gaussian average
+    // IntegrateGaussian does the actual
+    // integration of diffusion gaussian
+    IntegrateGaussian(k, par0[k], w, sigma[k], inf[k], sup[k],tav);
+    
+    //------------------------------------------------------
+    // end Pside
+    //-------------------------------------------------------
+    
+    //------------------------------------------------------
+    // Nside
+    //-------------------------------------------------------
+    k=1;
+
+    xp=x; zp=z; 
 
-    Float_t dx = (xE - xI);
-    Float_t dz = (zE - zI);
-              
+    // correction for the Lorentz's angle
+    if(fLorentz) {
+      Float_t deltaxn = ((seg->Dy()*1.0E-4)/2-y)*tanLorAngN;
+      xp+=deltaxn;
+    }
     
-    // Debuging
-    /*
-    fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
-   
-       printf("%5d %8.3f %8.3f %8.3f %8.3f %d %d  %d\n", 
-             hitNo, xI, zI, dx, dz, 
-             stripP, stripN, track);
-     printf("%10.5f %10d \n", ionization, hitI->fTrack); 
-    */ 
+
+    seg->GetPadTxz(xp,zp);
+
+    tdrift[1] = ((seg->Dy()*1.0E-4)/2-y)/GetDriftVelocity(1);
     
-    // end of debuging   
+    //tang[k]=TMath::Tan(tang[k]);
     
+    w = zp; // N side strip number
     
-    eP=0;
-    eN=0;
+    if((w<(-0.5)) || (w>(GetNStrips()-0.5))) {
+      // this check rejects hits in regions not covered by strips
+      // 0.5 takes into account boundaries 
+      if(GetDebug(4)) cout << "Dead SSD region, x,z="<<x<<","<<z<<endl;
+      return; // There are dead region on the SSD sensitive volume.
+    } // end if
     
-    for (i=0; i<fSteps; i++) {
-        
-      //        arrayEP[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
-      //        arrayEN[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
-        arrayEP[i] = ionization/fSteps;
-        arrayEN[i] = ionization/fSteps;
-       
-        eP += arrayEP[i];
-        eN += arrayEN[i];
-    } 
-       
-    const Float_t kconv = 1.0e9 / 3.6;  // GeV -> e-hole pairs
-       
-    for(i=0; i<fSteps; i++) {
+      // sigma is the standard deviation of the diffusion gaussian
+    if(tdrift[k]<0) return;
     
-        arrayEP[i] = kconv * arrayEP[i] * (ionization / eP); 
-        arrayEN[i] = kconv * arrayEN[i] * (ionization / eN);        
-    }    
-        
-    dx /= fSteps;
-    dz /= fSteps;  
-
-    Float_t sigmaP, sigmaN; 
-    fResponse->SigmaSpread(sigmaP,sigmaN); 
-
-    //printf("SigmaP SigmaN %f %f\n",sigmaP, sigmaN);
-
-    Float_t noiseP, noiseN;
-    fResponse->GetNoiseParam(noiseP,noiseN);
-
-    //printf("NoiseP NoiseN %f %f\n",noiseP, noiseN);
-
-    for(i=0; i<fSteps; i++) {
-        
-        Int_t j;
-       
-        fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
-        //printf("hitNo %d i xI zI stripP stripN %d %f %f %d %d\n",hitNo,i,xI, zI, stripP, stripN);
-        dsP    = Get2Strip(1,stripP,xI, zI); // Between 0-1
-        dsN    = Get2Strip(0,stripN,xI, zI); // Between 0-1
-
-        sP = sigmaP * sqrt(300. * i / (fSteps));
-        sN = sigmaN * sqrt(300. * i /(fSteps-i));
-
-
-        sP = (i<2        && dsP>0.3 && dsP<0.7)? 20. : sP;  // square of (microns) 
-        sN = (i>fSteps-2 && dsN>0.3 && dsN<0.7)? 20. : sN;  // square of (microns) 
-
-        sP = (i==2 && dsP>0.4 && dsP<0.6)? 15. : sP;  // square of (microns) 
-        sN = (i==8 && dsN>0.4 && dsN<0.6)? 15. : sN;  // square of (microns)        
-        
-        for (j=-1; j<2; j++) {
-            if (stripP+j<0 || stripP+j>fNstrips) continue;
-            signal = arrayEP[i] * TMath::Abs( (F(j+0.5-dsP,sP)-F(j-0.5-dsP,sP)) );
-            //printf("SimSSD::HitsToDigits:%d arrayEP[%d]=%e signal=%e\n",j,i,arrayEP[i],signal);
-           if (signal > noiseP/fSteps) {
-               (*fP)[stripP+j] += signal;
-               dict = (fTracksP+stripP+j);   
-               (*dict).AddTrack(track);
-            } 
-       }  // end for j loop over neighboring strips
-
-        for (j=-1; j<2; j++) {
-            if (stripN+j<0 || stripN+j>fNstrips) continue;
-            signal = arrayEN[i] * TMath::Abs( (F(j+0.5-dsN,sN)-F(j-0.5-dsN,sN)) );
-            //printf("SimSSD::HitsToDigits:%d arrayEN[%d]=%e signal=%e\n",j,i,arrayEN[i],signal);
-            if (signal > noiseN/fSteps) {
-               (*fN)[stripN+j] += signal;
-               dict = (fTracksN+stripN+j);    //co to jest
-               (*dict).AddTrack(track);
-            } 
-        }  // end for j loop over neighboring strips
-                
-        xI += dx; 
-        zI += dz; 
-    }
+    sigma[k] = TMath::Sqrt(2*GetDiffConst(k)*tdrift[k]);
+    sigma[k] /= (GetStripPitch()*1.0E-4);  //units of Pitch
+    
+    if(sigma[k]==0.0) {        
+      Error("HitToDigit"," sigma[%d]=0",k);
+      exit(0);
+    } // end if
     
+    par0[k] = pairs;
+    // we integrate the diffusion gaussian from -3sigma to 3sigma 
+    inf[k] = w - 3*sigma[k]; // 3 sigma from the gaussian average  
+    sup[k] = w + 3*sigma[k]; // 3 sigma from the gaussian average
+    // IntegrateGaussian does the actual
+    // integration of diffusion gaussian
+    IntegrateGaussian(k, par0[k], w, sigma[k], inf[k], sup[k],tav);
     
+    //-------------------------------------------------
+    // end Nside
+    //-------------------------------------------------
+    
+    
+  } // end stepping
 }
 
-
-//____________________________________________________________________
-//
-//  Private Methods for Simulation
 //______________________________________________________________________
-//
-
-void AliITSsimulationSSD::ApplyNoise() {
+void AliITSsimulationSSD::ApplyNoise(AliITSpList *pList,Int_t module){
   // Apply Noise.
-   Float_t noiseP, noiseN;
-   fResponse->GetNoiseParam(noiseP,noiseN);
-       
-    Int_t i;
-    for(i = 0; i<fNstrips; i++) {
-       (*fP)[i] += gRandom->Gaus(0,noiseP);
-       (*fN)[i] += gRandom->Gaus(0,noiseN);
-    }
+  Int_t ix;
+  Double_t signal,noise;
+  AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module);
+   
+  // Pside
+  for(ix=0;ix<GetNStrips();ix++){      // loop over strips
+    
+    // noise is gaussian
+    noise  = (Double_t) gRandom->Gaus(0,res->GetNoiseP(ix));
+    
+    // need to calibrate noise 
+    // NOTE. noise from the calibration database comes uncalibrated, 
+    // it needs to be calibrated in order to be added
+    // to the signal. It will be decalibrated later on together with the noise    
+    noise *= (Double_t) res->GetGainP(ix); 
+    
+    // noise comes in ADC channels from the calibration database
+    // It needs to be converted back to electronVolts
+    noise /= res->GetSSDDEvToADC(1.);
+    
+    // Finally, noise is added to the signal
+    signal = noise + fMapA2->GetSignal(0,ix);//get signal from map
+    fMapA2->SetHit(0,ix,signal); // give back signal to map
+    if(signal>0.0) pList->AddNoise(0,ix,module,noise);
+  } // loop over strip 
+  
+    // Nside
+  for(ix=0;ix<GetNStrips();ix++){      // loop over strips
+    noise  = (Double_t) gRandom->Gaus(0,res->GetNoiseN(ix));// give noise to signal
+    noise *= (Double_t) res->GetGainN(ix); 
+    noise /= res->GetSSDDEvToADC(1.);
+    signal = noise + fMapA2->GetSignal(1,ix);//get signal from map
+    fMapA2->SetHit(1,ix,signal); // give back signal to map
+    if(signal>0.0) pList->AddNoise(1,ix,module,noise);
+  } // loop over strip 
+  
 }
-
-//_________________________________________________________________________
-
-void AliITSsimulationSSD::ApplyCoupling() {
-  // Apply the effecto of electronic coupling between channels    
-    Int_t i;
-    for(i = 1; i<fNstrips-1; i++) {
-      (*fP)[i] += (*fP)[i-1]*fDCS->GetCouplingPL() + (*fP)[i+1]*fDCS->GetCouplingPR();
-      (*fN)[i] += (*fN)[i-1]*fDCS->GetCouplingNL() + (*fN)[i+1]*fDCS->GetCouplingNR();
-    }
+//______________________________________________________________________
+void AliITSsimulationSSD::ApplyCoupling(AliITSpList *pList,Int_t module) {
+  // Apply the effect of electronic coupling between channels
+  Int_t ix;
+  Double_t signal=0;
+  //AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module);
+  AliITSSimuParam* res = fDetType->GetSimuParam();
+    
+  Double_t *contrLeft  = new Double_t[GetNStrips()];
+  Double_t *contrRight = new Double_t[GetNStrips()];
+  
+  // P side coupling
+  for(ix=0;ix<GetNStrips();ix++){
+    if(ix>0) contrLeft[ix] = fMapA2->GetSignal(0,ix-1)*res->GetSSDCouplingPL();
+    else contrLeft[ix] = 0.0;
+    if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(0,ix+1)*res->GetSSDCouplingPR();
+    else contrRight[ix] = 0.0;
+  } // loop over strips 
+  
+  for(ix=0;ix<GetNStrips();ix++){
+    signal = contrLeft[ix] + contrRight[ix] - res->GetSSDCouplingPL() * fMapA2->GetSignal(0,ix)
+      - res->GetSSDCouplingPR() * fMapA2->GetSignal(0,ix);
+    fMapA2->AddSignal(0,ix,signal);
+    if(signal>0.0) pList->AddNoise(0,ix,module,signal);
+  } // loop over strips 
+  
+  // N side coupling
+  for(ix=0;ix<GetNStrips();ix++){
+    if(ix>0) contrLeft[ix] = fMapA2->GetSignal(1,ix-1)*res->GetSSDCouplingNL();
+    else contrLeft[ix] = 0.0;
+    if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(1,ix+1)*res->GetSSDCouplingNR();
+    else contrRight[ix] = 0.0;
+  } // loop over strips 
+  
+  for(ix=0;ix<GetNStrips();ix++){
+    signal = contrLeft[ix] + contrRight[ix] - res->GetSSDCouplingNL() * fMapA2->GetSignal(0,ix)
+      - res->GetSSDCouplingNR() * fMapA2->GetSignal(0,ix);
+    fMapA2->AddSignal(1,ix,signal);
+    if(signal>0.0) pList->AddNoise(1,ix,module,signal);
+  } // loop over strips 
+  
+
+  delete [] contrLeft;
+  delete [] contrRight; 
 }
 
-//__________________________________________________________________________
-
-void AliITSsimulationSSD::ApplyThreshold() {
-  // Applies the effect of a threshold on the signals for digitization.
-   Float_t noiseP, noiseN;
-   fResponse->GetNoiseParam(noiseP,noiseN);
-
-   // or introduce the SetThresholds in fResponse  
+//______________________________________________________________________
+void AliITSsimulationSSD::ApplyDeadChannels(Int_t module) {
+  // Kill dead channels setting gain to zero
 
-    Int_t i;
-    for(i=0; i<fNstrips; i++) {
-       (*fP)[i] = ((*fP)[i] > noiseP*4) ? (*fP)[i] : 0;
-       (*fN)[i] = ((*fN)[i] > noiseN*4) ? (*fN)[i] : 0; 
-    }
-        
-}
+  AliITSCalibrationSSD* res = (AliITSCalibrationSSD*)GetCalibrationModel(module);
 
-//__________________________________________________________________________
+  for(Int_t i=0;i<GetNStrips();i++){
 
-void AliITSsimulationSSD::ApplyDAQ() {
-  // Converts simulated signals to simulated ADC counts
-    AliITS *its=(AliITS*)gAlice->GetModule("ITS");
+    if(res->IsPChannelBad(i)) res->SetGainP(i,0.0);
+    if(res->IsNChannelBad(i)) res->SetGainN(i,0.0);
 
-    Float_t noiseP, noiseN;
-    fResponse->GetNoiseParam(noiseP,noiseN);
+  } // loop over strips 
 
-    char opt[30],dummy[20];
-    fResponse->ParamOptions(opt,dummy);
+}
 
-    Int_t i,j;
-    if (strstr(opt,"SetInvalid")) {
-      // Set signal = 0 if invalid strip
-      for(i=0; i<fNstrips; i++) {
-         if (!(fDCS->IsValidP(i))) (*fP)[i] = 0;
-        if (!(fDCS->IsValidN(i))) (*fN)[i] = 0;
-      }
-    }
-    
-    Int_t digits[3], tracks[3], hits[3];
-    Float_t charges[3];
-    Float_t phys=0;
-    for(i=0;i<3;i++) tracks[i]=0;
-    for(i=0; i<fNstrips; i++) { 
-       if( (strstr(opt,"SetInvalid") && (*fP)[i] < noiseP*4) || !(*fP)[i]) continue;
-          digits[0]=1;
-          digits[1]=i;
-          digits[2]=(int)(*fP)[i];
-          for(j=0; j<(fTracksP+i)->GetNTracks(); j++) {
-            if(j>2) continue;
-           if((fTracksP+i)->GetNTracks()) tracks[j]=(fTracksP+i)->GetTrack(j);
-           else tracks[j]=-2;
-           charges[j] = 0;
-           hits[j] = -1;
-          }
-          its->AddSimDigit(2,phys,digits,tracks,hits,charges);
-          
-    }
-    
-    
-    for(i=0; i<fNstrips; i++) {
-       if( (strstr(opt,"SetInvalid") && (*fN)[i] < noiseN*4)|| !(*fN)[i]) continue;
-          digits[0]=0;
-          digits[1]=i;
-          digits[2]=(int)(*fN)[i];
-          for( j=0; j<(fTracksN+i)->GetNTracks(); j++) {
-            if(j>2) continue;
-            if((fTracksN+i)->GetNTracks()) tracks[j]=(fTracksN+i)->GetTrack(j);
-           else tracks[j]=-2;
-            charges[j] = 0;
-           hits[j] = -1;
-          }
-          its->AddSimDigit(2,phys,digits,tracks,hits,charges);
-          
-    }
+//______________________________________________________________________
+Float_t AliITSsimulationSSD::F(Float_t av, Float_t x, Float_t s) {
+    // Computes the integral of a gaussian using Error Function
+    Float_t sqrt2 = TMath::Sqrt(2.0);
+    Float_t sigm2 = sqrt2*s;
+    Float_t integral;
+
+    integral = 0.5 * AliMathBase::ErfFast( (x - av) / sigm2);
+    return integral;
+}
+//______________________________________________________________________
+void AliITSsimulationSSD::IntegrateGaussian(Int_t k,Double_t par, Double_t w,
+                                            Double_t sigma, 
+                                            Double_t inf, Double_t sup,
+                                            AliITSTableSSD *tav) {
+    // integrate the diffusion gaussian
+    // remind: inf and sup are w-3sigma and w+3sigma
+    //         we could define them here instead of passing them
+    //         this way we are free to introduce asimmetry
+
+    Double_t a=0.0, b=0.0;
+    Double_t dXCharge1 = 0.0, dXCharge2 = 0.0;
+    // dXCharge1 and 2 are the charge to two neighbouring strips
+    // Watch that we only involve at least two strips
+    // Numbers greater than 2 of strips in a cluster depend on
+    //  geometry of the track and delta rays, not charge diffusion!   
+
+    Double_t strip = TMath::Floor(w);         // closest strip on the left
+
+    if ( TMath::Abs((strip - w)) < 0.5) { 
+        // gaussian mean is closer to strip on the left
+        a = inf;                         // integration starting point
+        if((strip+0.5)<=sup) {
+            // this means that the tail of the gaussian goes beyond
+            // the middle point between strips ---> part of the signal
+            // is given to the strip on the right
+            b = strip + 0.5;               // integration stopping point
+            dXCharge1 = F( w, b, sigma) - F(w, a, sigma);
+            dXCharge2 = F( w, sup, sigma) - F(w ,b, sigma); 
+        }else { 
+            // this means that all the charge is given to the strip on the left
+            b = sup;
+            dXCharge1 = 0.9973;   // gaussian integral at 3 sigmas
+            dXCharge2 = 0.0;
+        } // end if
+        dXCharge1 = par * dXCharge1;// normalize by mean of number of carriers
+        dXCharge2 = par * dXCharge2;
+
+        // for the time being, signal is the charge
+        // in ChargeToSignal signal is converted in ADC channel
+        fMapA2->AddSignal(k,(Int_t)strip,dXCharge1);
+        tav->Add(k,(Int_t)strip);
+        if(((Int_t) strip) < (GetNStrips()-1)) {
+            // strip doesn't have to be the last (remind: last=GetNStrips()-1)
+            // otherwise part of the charge is lost
+            fMapA2->AddSignal(k,((Int_t)strip+1),dXCharge2);
+            tav->Add(k,((Int_t)(strip+1)));
+        } // end if
+    }else{
+        // gaussian mean is closer to strip on the right
+        strip++;     // move to strip on the rigth
+        b = sup;     // now you know where to stop integrating
+        if((strip-0.5)>=inf) { 
+            // tail of diffusion gaussian on the left goes left of
+            // middle point between strips
+            a = strip - 0.5;        // integration starting point
+            dXCharge1 = F(w, b, sigma) - F(w, a, sigma);
+            dXCharge2 = F(w, a, sigma) - F(w, inf, sigma);
+        }else {
+            a = inf;
+            dXCharge1 = 0.9973;   // gaussian integral at 3 sigmas
+            dXCharge2 = 0.0;
+        } // end if
+        dXCharge1 = par * dXCharge1;    // normalize by means of carriers
+        dXCharge2 = par * dXCharge2;
+        // for the time being, signal is the charge
+        // in ChargeToSignal signal is converted in ADC channel
+        fMapA2->AddSignal(k,(Int_t)strip,dXCharge1);
+        tav->Add(k,(Int_t)strip);
+        if(((Int_t) strip) > 0) {
+            // strip doesn't have to be the first
+            // otherwise part of the charge is lost
+            fMapA2->AddSignal(k,((Int_t)strip-1),dXCharge2);
+            tav->Add(k,((Int_t)(strip-1)));
+        } // end if
+    } // end if
+}
+//______________________________________________________________________
+Int_t AliITSsimulationSSD::NumOfSteps(Double_t x, Double_t y, Double_t z,
+                                      Double_t &dex,Double_t &dey,
+                                      Double_t &dez){
+    // number of steps
+    // it also returns steps for each coord
+    //AliITSsegmentationSSD *seg = new AliITSsegmentationSSD();
+
+    Double_t step = 25E-4;
+    //step = (Double_t) seg->GetStepSize();  // step size (cm)
+    Int_t numOfSteps = (Int_t) (TMath::Sqrt(x*x+y*y+z*z)/step); 
+
+    if (numOfSteps < 1) numOfSteps = 1;       // one step, at least
+    //numOfSteps=1;
+
+    // we could condition the stepping depending on the incident angle
+    // of the track
+    dex = x/numOfSteps;
+    dey = y/numOfSteps;
+    dez = z/numOfSteps;
     
+    return numOfSteps;
 }
+//----------------------------------------------------------------------
+void AliITSsimulationSSD::GetList(Int_t label,Int_t hit,Int_t mod,
+                                  AliITSpList *pList,AliITSTableSSD *tav) {
+    // loop over nonzero digits
+    Int_t ix,i;
+    Double_t signal=0.;
+
+    for(Int_t k=0; k<2; k++) {
+        ix=tav->Use(k);
+        while(ix>-1){
+            signal = fMapA2->GetSignal(k,ix);
+            if(signal==0.0) {
+                ix=tav->Use(k);
+                continue;
+            } // end if signal==0.0
+            // check the signal magnitude
+            for(i=0;i<pList->GetNSignals(k,ix);i++){
+                signal -= pList->GetTSignal(k,ix,i);
+            } // end for i
+            //  compare the new signal with already existing list
+            if(signal>0)pList->AddSignal(k,ix,label,hit,mod,signal);
+            ix=tav->Use(k);
+        } // end of loop on strips
+    } // end of loop on P/N side
+    tav->Clear();
+}
+//----------------------------------------------------------------------
+void AliITSsimulationSSD::ChargeToSignal(Int_t module,const AliITSpList *pList) {
+    // charge to signal
+    static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
+    Float_t threshold = 0.;
+    Int_t size = AliITSdigitSSD::GetNTracks();
+    Int_t * digits = new Int_t[size];
+    Int_t * tracks = new Int_t[size];
+    Int_t * hits = new Int_t[size];
+    Int_t j1;
+    Float_t charges[3] = {0.0,0.0,0.0};
+    Float_t signal;
+    AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module);
+    AliITSSimuParam* simpar = fDetType->GetSimuParam();
+
+    for(Int_t k=0;k<2;k++){         // both sides (0=Pside, 1=Nside)
+      for(Int_t ix=0;ix<GetNStrips();ix++){     // loop over strips
+
+       // if strip is dead -> gain=0
+       if( ((k==0)&&(res->GetGainP(ix)==0)) || ((k==1)&&(res->GetGainN(ix)==0))) continue;
+       
+       signal = fMapA2->GetSignal(k,ix);
+       // signal has to be uncalibrated
+       // In real life, gains are supposed to be calculated from calibration runs,
+       // stored in the calibration DB and used in the reconstruction
+       // (see AliITSClusterFinderSSD.cxx)
+       if(k==0) signal /= res->GetGainP(ix);
+       else signal /= res->GetGainN(ix);
+
+       // signal is converted in unit of ADC
+       signal = res->GetSSDDEvToADC(signal);
+       if(signal>4095.) signal = 4095.;//if exceeding, accumulate last one
+
+       // threshold for zero suppression is set on the basis of the noise
+       // A good value is 3*sigma_noise
+       if(k==0) threshold = res->GetNoiseP(ix);
+       else threshold = res->GetNoiseN(ix);
+
+       threshold *= simpar->GetSSDZSThreshold(); // threshold at 3 sigma noise
+
+       if(signal < threshold) continue;
+       //cout<<signal<<" "<<threshold<<endl;
+
+       digits[0] = k;
+       digits[1] = ix;
+       digits[2] = TMath::Nint(signal);
+       for(j1=0;j1<size;j1++)if(j1<pList->GetNEntries()){
+         // only three in digit.
+         tracks[j1]  = pList->GetTrack(k,ix,j1);
+         hits[j1]    = pList->GetHit(k,ix,j1);
+       }else{
+         tracks[j1]  = -3;
+         hits[j1]    = -1;
+       } // end for j1
+       // finally add digit
+       aliITS->AddSimDigit(2,0,digits,tracks,hits,charges);
+      } // end for ix
+    } // end for k
+    delete [] digits;
+    delete [] tracks;
+    delete [] hits;
+}
+//______________________________________________________________________
+void AliITSsimulationSSD::WriteSDigits(AliITSpList *pList){
+    // Fills the Summable digits Tree
+    Int_t i,ni,j,nj;
+    static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
+
+    pList->GetMaxMapIndex(ni,nj);
+    for(i=0;i<ni;i++)for(j=0;j<nj;j++){
+        if(pList->GetSignalOnly(i,j)>0.0){
+            aliITS->AddSumDigit(*(pList->GetpListItem(i,j)));
+            if(GetDebug(4)) cout << "pListSSD: "<<*(pList->GetpListItem(i,j))
+                                << endl;
+        } // end if
+    } // end for i,j
+  return;
+}
+//______________________________________________________________________
+void AliITSsimulationSSD::FillMapFrompList(AliITSpList *pList){
+    // Fills fMap2A from the pList of Summable digits
+    Int_t k,ix;
 
-
-//____________________________________________________________________________
-
-Float_t AliITSsimulationSSD::F(Float_t x, Float_t s) {
-  // Computes the integral of a gaussian at the mean valuse x with sigma s.
-
-    //printf("SDD:F(%e,%e)\n",x,s);
-    return 0.5*TMath::Erf(x * fPitch / s) ;
-} 
-
+    for(k=0;k<2;k++)for(ix=0;ix<GetNStrips();ix++) 
+        fMapA2->AddSignal(k,ix,pList->GetSignal(k,ix));
+    return;
+}
 //______________________________________________________________________
+void AliITSsimulationSSD::Print(ostream *os){
+    //Standard output format for this class
+
+    //AliITSsimulation::Print(os);
+    *os << fIonE <<",";
+    *os << fDifConst[0] <<","<< fDifConst[1] <<",";
+    *os << fDriftVel[0] <<","<< fDriftVel[1];
+    //*os <<","; fDCS->Print(os);
+    //*os <<","; fMapA2->Print(os);
+}
+//______________________________________________________________________
+void AliITSsimulationSSD::Read(istream *is){
+    // Standard output streaming function.
+
+    //AliITSsimulation::Read(is);
+    *is >> fIonE;
+    *is >> fDifConst[0] >> fDifConst[1];
+    *is >> fDriftVel[0] >> fDriftVel[1];
+    //fDCS->Read(is);
+    //fMapA2->Read(is);
+}
+//______________________________________________________________________
+ostream &operator<<(ostream &os,AliITSsimulationSSD &source){
+    // Standard output streaming function.
 
-Float_t AliITSsimulationSSD::Get2Strip(Int_t flag, Int_t iStrip, Float_t x, Float_t z){
-  // Returns the relative space between two strips.
+    source.Print(&os);
+    return os;
+}
+//______________________________________________________________________
+istream &operator>>(istream &os,AliITSsimulationSSD &source){
+    // Standard output streaming function.
 
-    // flag==1 for Pside, 0 for Nside
+    source.Read(&os);
+    return os;
+}
+//______________________________________________________________________
 
-    Float_t stereoP, stereoN;
-    fSegmentation->Angles(stereoP,stereoN);
-    
-    Float_t tanP=TMath::Tan(stereoP);
-    Float_t tanN=TMath::Tan(stereoN);
-    Float_t dx = fSegmentation->Dx();
-    Float_t dz = fSegmentation->Dz();
 
 
-     x += dx/2;
-     z += dz/2; 
-    
-     if (flag) return (x - z*tanP) / fPitch - iStrip;       // from 0 to 1
-     else  return (x - tanN*(dz - z)) / fPitch - iStrip;
-}
-//____________________________________________________________________________