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