ITS/UPGRADE/AliITSUSimulationPix.cxx

   1 /*
   2  questions to experts: why RemoveDeadPixels should be called before FrompListToDigits ?
   3
   4
   5 */
   6
   7 /**************************************************************************
   8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   9 *                                                                        *
  10 * Author: The ALICE Off-line Project.                                    *
  11 * Contributors are mentioned in the code where appropriate.              *
  12 *                                                                        *
  13 * Permission to use, copy, modify and distribute this software and its   *
  14 * documentation strictly for non-commercial purposes is hereby granted   *
  15 * without fee, provided that the above copyright notice appears in all   *
  16 * copies and that both the copyright notice and this permission notice   *
  17 * appear in the supporting documentation. The authors make no claims     *
  18 * about the suitability of this software for any purpose. It is          *
  19 * provided "as is" without express or implied warranty.                  *
  20 **************************************************************************/
  21
  22 #include <TGeoGlobalMagField.h>
  23 #include <TH1.h>
  24 #include <TString.h>
  25 #include "AliITSU.h"
  26 #include "AliITSUDigitPix.h"
  27 #include "AliITSUHit.h"
  28 #include "AliITSUModule.h"
  29 #include "AliITSUSensMap.h"
  30 #include "AliITSUCalibrationPix.h"
  31 #include "AliITSUSegmentationPix.h"
  32 #include "AliITSUSimulationPix.h"
  33 #include "AliLog.h"
  34 #include "AliRun.h"
  35 #include "AliMagF.h"
  36 #include "AliMathBase.h"
  37 #include "AliITSUSimuParam.h"
  38 #include "AliITSUSDigit.h"
  39
  40 using std::cout;
  41 using std::endl;
  42
  43 ClassImp(AliITSUSimulationPix)
  44 ////////////////////////////////////////////////////////////////////////
  45 //  Version: 1
  46 //  Modified by D. Elia, G.E. Bruno, H. Tydesjo
  47 //  Fast diffusion code by Bjorn S. Nilsen
  48 //  March-April 2006
  49 //  October     2007: GetCalibrationObjects() removed
  50 //
  51 //  Version: 0
  52 //  Written by Boris Batyunya
  53 //  December 20 1999
  54 //
  55 //  Adapted for pixels of ITS upgrade July 2012, ruben.shahoyan@cern.ch
  56 //
  57 //  AliITSUSimulationPix is to do the simulation of pixels
  58 //
  59 ////////////////////////////////////////////////////////////////////////
  60
  61 //______________________________________________________________________
  62 AliITSUSimulationPix::AliITSUSimulationPix()
  63 :  fTanLorAng(0)
  64   ,fStrobe(kTRUE)
  65   ,fStrobeLenght(4)
  66   ,fStrobePhase(-12.5e-9)
  67 {
  68    // Default constructor.
  69 }
  70
  71 //______________________________________________________________________
  72 AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap* map)
  73   :AliITSUSimulation(sim,map)
  74   ,fTanLorAng(0)
  75   ,fStrobe(kTRUE)
  76   ,fStrobeLenght(4)
  77   ,fStrobePhase(-12.5e-9)
  78 {
  79   // standard constructor
  80   Init();
  81 }
  82
  83 //______________________________________________________________________
  84 AliITSUSimulationPix::AliITSUSimulationPix(const AliITSUSimulationPix &s)
  85   :AliITSUSimulation(s)
  86   ,fTanLorAng(s.fTanLorAng)
  87   ,fStrobe(s.fStrobe)
  88   ,fStrobeLenght(s.fStrobeLenght)
  89   ,fStrobePhase(s.fStrobePhase)
  90 {
  91   //     Copy Constructor
  92 }
  93
  94
  95 //______________________________________________________________________
  96 AliITSUSimulationPix::~AliITSUSimulationPix()
  97 {
  98   // destructor
  99   // only the sens map is owned and it is deleted by ~AliITSUSimulation
 100 }
 101
 102 //______________________________________________________________________
 103 AliITSUSimulationPix& AliITSUSimulationPix::operator=(const AliITSUSimulationPix &s)
 104 {
 105   //    Assignment operator
 106   if (&s == this) return *this;
 107   AliITSUSimulation::operator=(s);
 108   fStrobe       = s.fStrobe;
 109   fStrobeLenght = s.fStrobeLenght;
 110   fStrobePhase  = s.fStrobePhase;
 111   return *this;
 112 }
 113
 114 //______________________________________________________________________
 115 void AliITSUSimulationPix::Init()
 116 {
 117   // Initilization
 118   if (fSimuParam->GetPixLorentzDrift()) SetTanLorAngle(fSimuParam->GetPixLorentzHoleWeight());
 119 }
 120
 121 //______________________________________________________________________
 122 Bool_t AliITSUSimulationPix::SetTanLorAngle(Double_t weightHole)
 123 {
 124   // This function set the Tangent of the Lorentz angle.
 125   // A weighted average is used for electrons and holes
 126   // Input: Double_t weightHole: wheight for hole: it should be in the range [0,1]
 127   // output: Bool_t : kTRUE in case of success
 128   //
 129   if (weightHole<0) {
 130     weightHole=0.;
 131     AliWarning("You have asked for negative Hole weight");
 132     AliWarning("I'm going to use only electrons");
 133   }
 134   if (weightHole>1) {
 135     weightHole=1.;
 136     AliWarning("You have asked for weight > 1");
 137     AliWarning("I'm going to use only holes");
 138   }
 139   Double_t weightEle=1.-weightHole;
 140   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 141   if (!fld) AliFatal("The field is not initialized");
 142   Double_t bz = fld->SolenoidField();
 143   fTanLorAng = TMath::Tan(weightHole*fSimuParam->LorentzAngleHole(bz) +
 144                           weightEle*fSimuParam->LorentzAngleElectron(bz));
 145   return kTRUE;
 146 }
 147
 148 //_____________________________________________________________________
 149 void AliITSUSimulationPix::SDigitiseModule(AliITSUModule *mod, Int_t /*mask*/,Int_t event, AliITSsegmentation* seg)
 150 {
 151   //  This function begins the work of creating S-Digits.
 152   if (!(mod->GetNHits())) {
 153     AliDebug(1,Form("In event %d module %d there are %d hits returning.",
 154                     event, mod->GetIndex(),mod->GetNHits()));
 155     return;
 156   }
 157   //
 158   if (fStrobe) if (event != fEvent) GenerateStrobePhase();
 159   InitSimulationModule(mod->GetIndex(), event, seg );
 160   // Hits2SDigits(mod);
 161   Hits2SDigitsFast(mod);
 162   if (fSimuParam->GetPixAddNoisyFlag())   AddNoisyPixels();
 163   if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();
 164   WriteSDigits();
 165   ClearMap();
 166 }
 167
 168 //______________________________________________________________________
 169 void AliITSUSimulationPix::WriteSDigits()
 170 {
 171   //  This function adds each S-Digit to pList
 172   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
 173   int nsd = fSensMap->GetEntries();
 174   for (int i=0;i<nsd;i++) {
 175     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
 176     if (!sd->GetSumSignal()>0 || fSensMap->IsDisabled(sd)) continue;
 177     aliITS->AddSumDigit(*sd);
 178   }
 179   return;
 180 }
 181
 182 //______________________________________________________________________
 183 void AliITSUSimulationPix::FinishSDigitiseModule()
 184 {
 185    //  This function calls SDigitsToDigits which creates Digits from SDigits
 186   FrompListToDigits();
 187   ClearMap();
 188   return;
 189 }
 190
 191 //______________________________________________________________________
 192 void AliITSUSimulationPix::DigitiseModule(AliITSUModule *mod,Int_t /*mask*/, Int_t event, AliITSsegmentation* seg)
 193 {
 194   //  This function creates Digits straight from the hits and then adds
 195   //  electronic noise to the digits before adding them to pList
 196   //  Each of the input variables is passed along to Hits2SDigits
 197   //
 198   if (fStrobe) if (event != fEvent) GenerateStrobePhase();
 199   InitSimulationModule( mod->GetIndex(), event, seg );
 200   // Hits2SDigits(mod);
 201   Hits2SDigitsFast(mod);
 202   //
 203   if (fSimuParam->GetPixAddNoisyFlag())   AddNoisyPixels();
 204   if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();
 205   FrompListToDigits();
 206   ClearMap();
 207 }
 208
 209 //______________________________________________________________________
 210 void AliITSUSimulationPix::Hits2SDigits(AliITSUModule *mod)
 211 {
 212   // Does the charge distributions using Gaussian diffusion charge charing.
 213   const Double_t kBunchLenght = 25e-9; // LHC clock
 214   Int_t nhits = mod->GetNHits();
 215   if (!nhits) return;
 216   //
 217   Int_t h,ix,iz,i;
 218   Int_t idtrack;
 219   Float_t x,y,z;  // keep coordinates float (required by AliSegmentation)
 220   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;
 221   Double_t t,tp,st,dt=0.2,el,sig,sigx,sigz,fda;
 222   Double_t thick = 0.5*fSeg->Dy();  // Half Thickness
 223   fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
 224   //
 225   for (h=0;h<nhits;h++) {
 226     if (fStrobe &&
 227         ((mod->GetHit(h)->GetTOF()<fStrobePhase) ||
 228          (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
 229         ) continue;
 230     //
 231     if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
 232     st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
 233     if (st>0.0) {
 234       st = (Double_t)((Int_t)(st*1e4)); // number of microns
 235       if (st<=1.0) st = 1.0;
 236       dt = 1.0/st;               // RS TODO: do we need 1 micron steps?
 237       for (t=0.0;t<1.0;t+=dt) { // Integrate over t
 238         tp  = t+0.5*dt;
 239         x   = x0+x1*tp;
 240         y   = y0+y1*tp;
 241         z   = z0+z1*tp;
 242         if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 243         el  = dt * de / fSimuParam->GetGeVToCharge();
 244         //
 245         sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
 246         sigx=sig;
 247         sigz=sig*fda;
 248         if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
 249         SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
 250       } // end for t
 251     } else { // st == 0.0 deposit it at this point
 252       x   = x0;
 253       y   = y0;
 254       z   = z0;
 255       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 256       el  = de / fSimuParam->GetGeVToCharge();
 257       sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
 258       sigx = sig;
 259       sigz = sig*fda;
 260       if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
 261       SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
 262     } // end if st>0.0
 263   } // Loop over all hits h
 264   //
 265   // Coupling
 266   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
 267   AliITSUSDigit* dg = 0;
 268   switch (fSimuParam->GetPixCouplingOption()) {
 269   case AliITSUSimuParam::kNewCouplingPix :
 270     for (i=nd;i--;) {
 271       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 272       if (fSensMap->IsDisabled(dg)) continue;
 273       SetCoupling(dg,idtrack,h);
 274     }
 275     break;
 276   case AliITSUSimuParam::kOldCouplingPix:
 277     for (i=nd;i--;) {
 278       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 279       if (fSensMap->IsDisabled(dg)) continue;
 280       SetCouplingOld(dg,idtrack,h);
 281     }
 282     break;
 283   default:
 284     break;
 285
 286   } // end switch
 287   if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
 288 }
 289
 290 //______________________________________________________________________
 291 void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod)
 292 {
 293   // Does the charge distributions using Gaussian diffusion charge charing.    // Inputs:
 294   //    AliITSUModule *mod  Pointer to this module
 295   // Output:
 296   //    none.
 297   // Return:
 298   //    none.
 299   const Int_t kn10=10;
 300   const Double_t kti[kn10]={7.443716945e-3,2.166976971e-1,3.397047841e-1,
 301                             4.325316833e-1,4.869532643e-1,5.130467358e-1,
 302                             5.674683167e-1,6.602952159e-1,7.833023029e-1,
 303                             9.255628306e-1};
 304   const Double_t kwi[kn10]={1.477621124e-1,1.346333597e-1,1.095431813e-1,
 305                             7.472567455e-2,3.333567215e-2,3.333567215e-2,
 306                             7.472567455e-2,1.095431813e-1,1.346333597e-1,
 307                             1.477621124e-1};
 308   const Double_t kBunchLenght = 25e-9; // LHC clock
 309   TObjArray *hits = mod->GetHits();
 310   Int_t nhits = hits->GetEntriesFast();
 311   if (nhits<=0) return;
 312   //
 313   Int_t h,ix,iz,i;
 314   Int_t idtrack;
 315   Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
 316   Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0;
 317   Double_t t,st,el,sig,sigx,sigz,fda,de=0.0,ld=0.0;
 318   Double_t thick = 0.5*fSeg->Dy();  // Half thickness
 319   fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
 320   //
 321   for (h=0;h<nhits;h++) {
 322     // Check if the hit is inside readout window
 323     if (fStrobe &&
 324         ((mod->GetHit(h)->GetTOF()<fStrobePhase) ||
 325          (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
 326         ) continue;
 327     //
 328     if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
 329     st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
 330     if (st>0.0)
 331       for (i=0;i<kn10;i++) { // Integrate over t
 332         t   = kti[i];
 333         x   = x0+x1*t;
 334         y   = y0+y1*t;
 335         z   = z0+z1*t;
 336         if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 337         el  = kwi[i]*de/fSimuParam->GetGeVToCharge();
 338         sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
 339         sigx=sig;
 340         sigz=sig*fda;
 341         if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
 342         SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
 343       } // end for i // End Integrate over t
 344     else { // st == 0.0 deposit it at this point
 345       x   = x0;
 346       y   = y0;
 347       z   = z0;
 348       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 349       el  = de / fSimuParam->GetGeVToCharge();
 350       sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
 351       sigx=sig;
 352       sigz=sig*fda;
 353       if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
 354       SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
 355     } // end if st>0.0
 356
 357   } // Loop over all hits h
 358
 359   // Coupling
 360   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
 361   AliITSUSDigit* dg = 0;
 362   switch (fSimuParam->GetPixCouplingOption()) {
 363   case AliITSUSimuParam::kNewCouplingPix :
 364     for (i=nd;i--;) {
 365       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 366       if (fSensMap->IsDisabled(dg)) continue;
 367       SetCoupling(dg,idtrack,h);
 368     }
 369   case AliITSUSimuParam::kOldCouplingPix:
 370     for (i=nd;i--;) {
 371       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 372       if (fSensMap->IsDisabled(dg)) continue;
 373       SetCouplingOld(dg,idtrack,h);
 374     }
 375     break;
 376   default:
 377     break;
 378   } // end switch
 379   if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
 380 }
 381
 382 //______________________________________________________________________
 383 void AliITSUSimulationPix::SpreadCharge(Double_t x0,Double_t z0,
 384                                           Int_t ix0,Int_t iz0,
 385                                           Double_t el,Double_t sig,Double_t ld,
 386                                           Int_t t,Int_t hi)
 387 {
 388    // Spreads the charge over neighboring cells. Assume charge is distributed
 389    // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg)
 390    // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig)
 391    // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
 392    // Defined this way, the integral over all x and z is el.
 393    // Inputs:
 394    //    Double_t x0   x position of point where charge is liberated
 395    //    Double_t z0   z position of point where charge is liberated
 396    //    Int_t    ix0  row of cell corresponding to point x0
 397    //    Int_t    iz0  columb of cell corresponding to point z0
 398    //    Double_t el   number of electrons liberated in this step
 399    //    Double_t sig  Sigma difusion for this step (y0 dependent)
 400    //    Double_t ld   lorentz drift in x for this step (y0 dependent)
 401    //    Int_t    t    track number
 402    //    Int_t    ti   hit track index number
 403    //    Int_t    hi   hit "hit" index number
 404    // Outputs:
 405    //     none.
 406    // Return:
 407    //     none.
 408    const Int_t knx = 3,knz = 2;
 409    const Double_t kRoot2 = 1.414213562; // Sqrt(2).
 410    Int_t ix,iz,ixs,ixe,izs,ize;
 411    Float_t x,z;  // keep coordinates float (required by AliSegmentation)
 412    Double_t s,sp,x1,x2,z1,z2;
 413    //
 414    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));
 415    if (sig<=0.0) { // if sig<=0 No diffusion to simulate.
 416      UpdateMapSignal(iz0,ix0,t,hi,el);
 417      return;
 418    } // end if
 419    sp = 1.0/(sig*kRoot2);
 420    ixs = TMath::Max(-knx+ix0,0);
 421    ixe = TMath::Min(knx+ix0,fSeg->Npx()-1);
 422    izs = TMath::Max(-knz+iz0,0);
 423    ize = TMath::Min(knz+iz0,fSeg->Npz()-1);
 424    for (ix=ixs;ix<=ixe;ix++)
 425      for (iz=izs;iz<=ize;iz++) {
 426        fSeg->DetToLocal(ix,iz,x,z); // pixel center
 427        double dxi = 0.5*fSeg->Dpx(ix);
 428        double dzi = 0.5*fSeg->Dpz(iz);
 429        x1  = x;
 430        z1  = z;
 431        x2  = x1 + dxi; // Upper
 432        x1 -= dxi;  // Lower
 433        z2  = z1 + dzi; // Upper
 434        z1 -= dzi;  // Lower
 435        x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
 436        x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
 437        z1 -= z0; // Distance from where track traveled
 438        z2 -= z0; // Distance from where track traveled
 439        s   = el*0.25; // Correction based on definision of Erfc
 440        s  *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(sp*x2);
 441        s  *= AliMathBase::ErfcFast(sp*z1) - AliMathBase::ErfcFast(sp*z2);
 442        if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
 443      } // end for ix, iz
 444    //
 445 }
 446
 447 //______________________________________________________________________
 448 void AliITSUSimulationPix::SpreadChargeAsym(Double_t x0,Double_t z0,
 449                                             Int_t ix0,Int_t iz0,
 450                                             Double_t el,Double_t sigx,Double_t sigz,
 451                                             Double_t ld,Int_t t,Int_t hi)
 452 {
 453   // Spreads the charge over neighboring cells. Assume charge is distributed
 454   // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
 455   // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
 456   // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
 457   // Defined this way, the integral over all x and z is el.
 458   // Inputs:
 459   //    Double_t x0   x position of point where charge is liberated
 460   //    Double_t z0   z position of point where charge is liberated
 461   //    Int_t    ix0  row of cell corresponding to point x0
 462   //    Int_t    iz0  columb of cell corresponding to point z0
 463   //    Double_t el   number of electrons liberated in this step
 464   //    Double_t sigx Sigma difusion along x for this step (y0 dependent)
 465   //    Double_t sigz Sigma difusion along z for this step (z0 dependent)
 466   //    Double_t ld   lorentz drift in x for this stip (y0 dependent)
 467   //    Int_t    t    track number
 468   //    Int_t    ti   hit track index number
 469   //    Int_t    hi   hit "hit" index number
 470   // Outputs:
 471   //     none.
 472   // Return:
 473   //     none.
 474   const Int_t knx = 3,knz = 3; // RS: TO TUNE
 475   const Double_t kRoot2 = 1.414213562; // Sqrt(2).
 476   Int_t ix,iz,ixs,ixe,izs,ize;
 477   Float_t x,z;   // keep coordinates float (required by AliSegmentation)
 478   Double_t s,spx,spz,x1,x2,z1,z2;
 479   //
 480   if (GetDebug(2)) AliInfo(Form("(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,sigx=%e, sigz=%e, t=%d,i=%d,ld=%e)",
 481                                 x0,z0,ix0,iz0,el,sigx,sigz,t,hi,ld));
 482   if (sigx<=0.0 || sigz<=0.0) { // if sig<=0 No diffusion to simulate.
 483     UpdateMapSignal(iz0,ix0,t,hi,el);
 484     return;
 485   } // end if
 486   spx = 1.0/(sigx*kRoot2);
 487   spz = 1.0/(sigz*kRoot2);
 488   ixs = TMath::Max(-knx+ix0,0);
 489   ixe = TMath::Min(knx+ix0,fSeg->Npx()-1);
 490   izs = TMath::Max(-knz+iz0,0);
 491   ize = TMath::Min(knz+iz0,fSeg->Npz()-1);
 492   for (ix=ixs;ix<=ixe;ix++)
 493     for (iz=izs;iz<=ize;iz++) {
 494       fSeg->DetToLocal(ix,iz,x,z); // pixel center
 495       double dxi = 0.5*fSeg->Dpx(ix);
 496       double dzi = 0.5*fSeg->Dpz(iz);
 497       x1  = x;
 498       z1  = z;
 499       x2  = x1 + dxi; // Upper
 500       x1 -= dxi;      // Lower
 501       z2  = z1 + dzi; // Upper
 502       z1 -= dzi;      // Lower
 503       x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
 504       x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
 505       z1 -= z0; // Distance from where track traveled
 506       z2 -= z0; // Distance from where track traveled
 507       s   = el*0.25; // Correction based on definision of Erfc
 508       s  *= AliMathBase::ErfcFast(spx*x1) - AliMathBase::ErfcFast(spx*x2);
 509       s  *= AliMathBase::ErfcFast(spz*z1) - AliMathBase::ErfcFast(spz*z2);
 510       if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
 511     } // end for ix, iz
 512   //
 513 }
 514
 515 //______________________________________________________________________
 516 void AliITSUSimulationPix::RemoveDeadPixels()
 517 {
 518   // Removes dead pixels on each module (ladder)
 519   // This should be called before going from sdigits to digits (FrompListToDigits)
 520
 521   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibDead();
 522   if (!calObj) return;
 523   //
 524   if (calObj->IsBad()) {ClearMap(); return;} // whole module is masked
 525   //
 526   // remove single bad pixels one by one
 527   int nsingle = calObj->GetNrBadSingle();
 528   UInt_t col,row;
 529   for (int i=nsingle;i--;) {
 530     calObj->GetBadPixelSingle(i,row,col);
 531     fSensMap->DeleteItem(col,row);
 532   }
 533   int nsd = fSensMap->GetEntriesUnsorted();
 534   for (int isd=nsd;isd--;) {
 535     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->AtUnsorted(isd);
 536     if (fSensMap->IsDisabled(sd)) continue;
 537     fSensMap->GetMapIndex(sd->GetUniqueID(),col,row);
 538     int chip = fSeg->GetChipFromChannel(0,col);
 539     //    if (calObj->IsChipMarkedBad(chip)) fSensMap->Disable(sd); // this will simple mark the hit as bad
 540     if (calObj->IsChipMarkedBad(chip)) fSensMap->DeleteItem(sd); // this will suppress hit in the sorted list
 541   }
 542   //
 543 }
 544
 545 //______________________________________________________________________
 546 void AliITSUSimulationPix::AddNoisyPixels()
 547 {
 548   // Adds noisy pixels on each module (ladder)
 549   // This should be called before going from sdigits to digits (FrompListToDigits)
 550   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibNoisy();
 551   if (!calObj) return;
 552   for (Int_t i=calObj->GetNrBad(); i--;) UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i),
 553                                                         10*fSimuParam->GetPixThreshold(fModule));
 554   //
 555 }
 556
 557 //______________________________________________________________________
 558 void AliITSUSimulationPix::FrompListToDigits()
 559 {
 560   // add noise and electronics, perform the zero suppression and add the
 561   // digit to the list
 562   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
 563   UInt_t ix,iz;
 564   Double_t sig;
 565   const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
 566   static AliITSUDigitPix dig;
 567   // RS: in principle:
 568   // 1) for every pixel w/o hit we have to generate a noise and activate the pixel if the noise exceeds the threshold.
 569   // 2) for every pixel with hit we should add random noise and check if the total signal exceeds the threshold
 570   // With many channels this will be too time consuming, hence I do the following
 571   // 1) Precalculate the probability that the nois alone will exceed the threshold.
 572   // 2) Chose randomly empty pixels according to this probability and apply the noise above threshold.
 573   // 3) For pixels having a hits apply the usual noise and compare with threshold
 574   //
 575   // RS may use for ordered random sample generation dl.acm.org/ft_gateway.cfm?id=356313&type=pdf
 576   //
 577   int maxInd = fSensMap->GetMaxIndex();
 578   double minProb = 0.1/maxInd;
 579   //
 580   int nsd = fSensMap->GetEntries();
 581   Int_t prevID=0,curID=0;
 582   TArrayI ordSampleInd(100),ordSample(100);
 583   //
 584   double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(fModule);
 585   fSimuParam->GetPixNoise(fModule, noiseSig, noiseMean);
 586   probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold
 587   //
 588   for (int i=0;i<nsd;i++) {
 589     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
 590     if (fSensMap->IsDisabled(sd)) continue;
 591     curID = (int)sd->GetUniqueID();
 592     //
 593     if (probNoisy>minProb) { // generate randomly noisy pixels above the threshold, with ID's between previous hit and current
 594       CreateNoisyDigits(prevID,curID,probNoisy, noiseSig, noiseMean);
 595       prevID = curID+1;
 596     }
 597     //
 598     if ((sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(fModule)) continue;
 599     if (TMath::Abs(sig)>2147483647.0) { //RS?
 600       //PH 2147483647 is the max. integer
 601       //PH This apparently is a problem which needs investigation
 602       AliWarning(Form("Too big or too small signal value %f",sig));
 603     }
 604     fSensMap->GetMapIndex(sd->GetUniqueID(),iz,ix);
 605     dig.SetCoord1(iz);
 606     dig.SetCoord2(ix);
 607     dig.SetSignal(1);
 608     dig.SetSignalPix((Int_t)sig);
 609     int ntr = sd->GetNTracks();
 610     for (int j=0;j<ntr;j++) {
 611       dig.SetTrack(j,sd->GetTrack(j));
 612       dig.SetHit(j,sd->GetHit(j));
 613     }
 614     for (int j=ntr;j<knmaxtrk;j++) {
 615       dig.SetTrack(j,-3);
 616       dig.SetHit(j,-1);
 617     }
 618     aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix, &dig);
 619   }
 620   // if needed, add noisy pixels with id from last real hit to maxID
 621   if (probNoisy>minProb) CreateNoisyDigits(prevID,maxInd,probNoisy, noiseSig, noiseMean);
 622   //
 623 }
 624
 625 //______________________________________________________________________
 626 Int_t AliITSUSimulationPix::CreateNoisyDigits(Int_t minID,Int_t maxID,double probNoisy, double noise, double base)
 627 {
 628   // create random noisy digits above threshold within id range [minID,maxID[
 629   // see FrompListToDigits for details
 630   //
 631   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
 632   UInt_t ix,iz;
 633   static AliITSUDigitPix dig;
 634   static TArrayI ordSampleInd(100),ordSample(100); //RS!!! static is not thread-safe!!!
 635   const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
 636   //
 637   Int_t ncand = 0;
 638   int npix = maxID-minID;
 639   if (npix<1 || (ncand=gRandom->Poisson(npix*probNoisy))<1) return ncand; // decide how many noisy pixels will be added
 640   ncand = GenOrderedSample(npix,ncand,ordSample,ordSampleInd);
 641   int* ordV = ordSample.GetArray();
 642   int* ordI = ordSampleInd.GetArray();
 643   for (int j=0;j<ncand;j++) {
 644     fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],iz,ix);   // create noisy digit
 645     dig.SetCoord1(iz);
 646     dig.SetCoord2(ix);
 647     dig.SetSignal(1);
 648     dig.SetSignalPix((Int_t)AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,base,noise));
 649     for (int k=knmaxtrk;k--;) {
 650       dig.SetTrack(k,-3);
 651       dig.SetHit(k,-1);
 652     }
 653     aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix,&dig);
 654     if (GetDebug(2)) AliInfo(Form("Add noisy pixel %d(%d/%d) Noise=%d",ordV[ordI[j]],iz,ix,dig.GetSignalPix()));
 655   }
 656   return ncand;
 657 }
 658
 659 //______________________________________________________________________
 660 void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old, Int_t ntrack, Int_t idhit)
 661 {
 662   //  Take into account the coupling between adiacent pixels.
 663   //  The parameters probcol and probrow are the probability of the
 664   //  signal in one pixel shared in the two adjacent pixels along
 665   //  the column and row direction, respectively.
 666   //  Note pList is goten via GetMap() and module is not need any more.
 667   //  Otherwise it is identical to that coded by Tiziano Virgili (BSN).
 668   //Begin_Html
 669   /*
 670     <img src="picts/ITS/barimodel_3.gif">
 671      </pre>
 672      <br clear=left>
 673      <font size=+2 color=red>
 674      <a href="mailto:tiziano.virgili@cern.ch"></a>.
 675      </font>
 676      <pre>
 677    */
 678    //End_Html
 679    // Inputs:
 680   // old                  existing AliITSUSDigit
 681   // Int_t ntrack         track incex number
 682   // Int_t idhit          hit index number
 683   UInt_t col,row;
 684   Double_t pulse1,pulse2;
 685   Double_t couplR=0.0,couplC=0.0;
 686   Double_t xr=0.;
 687   //
 688   fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
 689   fSimuParam->GetPixCouplingParam(couplC,couplR);
 690   if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
 691                                 col,row,ntrack,idhit,couplC,couplR));
 692   pulse2 = pulse1 = old->GetSignal();
 693   if (pulse1<fSimuParam->GetPixMinElToAdd()) return; // too small signal
 694   for (Int_t isign=-1;isign<=1;isign+=2) {
 695     //
 696     // loop in col direction
 697     int j1 = int(col) + isign;
 698     xr = gRandom->Rndm();
 699     if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
 700     //
 701     // loop in row direction
 702     int j2 = int(row) + isign;
 703     xr = gRandom->Rndm();
 704     if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
 705   }
 706   //
 707 }
 708
 709 //______________________________________________________________________
 710 void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, Int_t ntrack,Int_t idhit)
 711 {
 712   //  Take into account the coupling between adiacent pixels.
 713   //  The parameters probcol and probrow are the fractions of the
 714   //  signal in one pixel shared in the two adjacent pixels along
 715   //  the column and row direction, respectively.
 716   //Begin_Html
 717   /*
 718     <img src="picts/ITS/barimodel_3.gif">
 719     </pre>
 720     <br clear=left>
 721     <font size=+2 color=red>
 722     <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
 723     </font>
 724     <pre>
 725   */
 726   //End_Html
 727   // Inputs:
 728   // old            existing AliITSUSDigit
 729   // ntrack         track incex number
 730   // idhit          hit index number
 731   // module         module number
 732   //
 733   UInt_t col,row;
 734   Double_t pulse1,pulse2;
 735   Double_t couplR=0.0,couplC=0.0;
 736   //
 737   fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
 738   fSimuParam->GetPixCouplingParam(couplC,couplR);
 739   if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
 740                                 col,row,ntrack,idhit,couplC,couplR));
 741  //
 742  if (old->GetSignal()<fSimuParam->GetPixMinElToAdd()) return; // too small signal
 743  for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
 744    pulse2 = pulse1 = old->GetSignal();
 745    //
 746    int j1 = int(col)+isign;
 747    pulse1 *= couplC;
 748    if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(fModule))) pulse1 = old->GetSignal();
 749    else UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
 750
 751    // loop in row direction
 752    int j2 = int(row) + isign;
 753    pulse2 *= couplR;
 754    if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(fModule))) pulse2 = old->GetSignal();
 755    else UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
 756  } // for isign
 757 }
 758
 759 //______________________________________________________________________
 760 void AliITSUSimulationPix::GenerateStrobePhase()
 761 {
 762   // Generate randomly the strobe
 763   // phase w.r.t to the LHC clock
 764   // Done once per event
 765   const Double_t kBunchLenght = 25e-9; // LHC clock
 766   fStrobePhase = ((Double_t)gRandom->Integer(fStrobeLenght))*kBunchLenght-
 767     (Double_t)fStrobeLenght*kBunchLenght+
 768     kBunchLenght/2;
 769 }
 770