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 kmictocm = 1.0e-4; // convert microns to cm.
 214   const Double_t kcmtomic = 1.e4;
 215   const Double_t kBunchLenght = 25e-9; // LHC clock
 216   Int_t nhits = mod->GetNHits();
 217   if (!nhits) return;
 218   //
 219   Int_t h,ix,iz,i;
 220   Int_t idtrack;
 221   Float_t x,y,z;  // keep coordinates float (required by AliSegmentation)
 222   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;
 223   Double_t t,tp,st,dt=0.2,el,sig,sigx,sigz,fda;
 224   Double_t thick = 0.5*kmictocm*fSeg->Dy();  // Half Thickness
 225   fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
 226   //
 227   for (h=0;h<nhits;h++) {
 228     if (fStrobe &&
 229         ((mod->GetHit(h)->GetTOF()<fStrobePhase) ||
 230          (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
 231         ) continue;
 232     //
 233     if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
 234     st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
 235     if (st>0.0) {
 236       st = (Double_t)((Int_t)(st*kcmtomic)); // number of microns
 237       if (st<=1.0) st = 1.0;
 238       dt = 1.0/st;
 239       for (t=0.0;t<1.0;t+=dt) { // Integrate over t
 240         tp  = t+0.5*dt;
 241         x   = x0+x1*tp;
 242         y   = y0+y1*tp;
 243         z   = z0+z1*tp;
 244         if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 245         el  = dt * de / fSimuParam->GetGeVToCharge();
 246         //
 247         sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
 248         sigx=sig;
 249         sigz=sig*fda;
 250         if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
 251         SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
 252       } // end for t
 253     } else { // st == 0.0 deposit it at this point
 254       x   = x0;
 255       y   = y0;
 256       z   = z0;
 257       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 258       el  = de / fSimuParam->GetGeVToCharge();
 259       sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
 260       sigx = sig;
 261       sigz = sig*fda;
 262       if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
 263       SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
 264     } // end if st>0.0
 265   } // Loop over all hits h
 266   //
 267   // Coupling
 268   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
 269   AliITSUSDigit* dg = 0;
 270   switch (fSimuParam->GetPixCouplingOption()) {
 271   case AliITSUSimuParam::kNewCouplingPix :
 272     for (i=nd;i--;) {
 273       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 274       if (fSensMap->IsDisabled(dg)) continue;
 275       SetCoupling(dg,idtrack,h);
 276     }
 277     break;
 278   case AliITSUSimuParam::kOldCouplingPix:
 279     for (i=nd;i--;) {
 280       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 281       if (fSensMap->IsDisabled(dg)) continue;
 282       SetCouplingOld(dg,idtrack,h);
 283     }
 284     break;
 285   default:
 286     break;
 287
 288   } // end switch
 289   if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
 290 }
 291
 292 //______________________________________________________________________
 293 void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod)
 294 {
 295   // Does the charge distributions using Gaussian diffusion charge charing.    // Inputs:
 296   //    AliITSUModule *mod  Pointer to this module
 297   // Output:
 298   //    none.
 299   // Return:
 300   //    none.
 301   const Double_t kmictocm = 1.0e-4; // convert microns to cm.
 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   const Double_t kBunchLenght = 25e-9; // LHC clock
 312   TObjArray *hits = mod->GetHits();
 313   Int_t nhits = hits->GetEntriesFast();
 314   if (nhits<=0) return;
 315   //
 316   Int_t h,ix,iz,i;
 317   Int_t idtrack;
 318   Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
 319   Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0;
 320   Double_t t,st,el,sig,sigx,sigz,fda,de=0.0,ld=0.0;
 321   Double_t thick = 0.5*kmictocm*fSeg->Dy();  // Half thickness
 322   fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
 323   //
 324   for (h=0;h<nhits;h++) {
 325     // Check if the hit is inside readout window
 326     if (fStrobe &&
 327         ((mod->GetHit(h)->GetTOF()<fStrobePhase) ||
 328          (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
 329         ) continue;
 330     //
 331     if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
 332     st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
 333     if (st>0.0)
 334       for (i=0;i<kn10;i++) { // Integrate over t
 335         t   = kti[i];
 336         x   = x0+x1*t;
 337         y   = y0+y1*t;
 338         z   = z0+z1*t;
 339         if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 340         el  = kwi[i]*de/fSimuParam->GetGeVToCharge();
 341         sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
 342         sigx=sig;
 343         sigz=sig*fda;
 344         if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
 345         SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
 346       } // end for i // End Integrate over t
 347     else { // st == 0.0 deposit it at this point
 348       x   = x0;
 349       y   = y0;
 350       z   = z0;
 351       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 352       el  = de / fSimuParam->GetGeVToCharge();
 353       sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
 354       sigx=sig;
 355       sigz=sig*fda;
 356       if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
 357       SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
 358     } // end if st>0.0
 359
 360   } // Loop over all hits h
 361
 362   // Coupling
 363   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
 364   AliITSUSDigit* dg = 0;
 365   switch (fSimuParam->GetPixCouplingOption()) {
 366   case AliITSUSimuParam::kNewCouplingPix :
 367     for (i=nd;i--;) {
 368       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 369       if (fSensMap->IsDisabled(dg)) continue;
 370       SetCoupling(dg,idtrack,h);
 371     }
 372   case AliITSUSimuParam::kOldCouplingPix:
 373     for (i=nd;i--;) {
 374       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 375       if (fSensMap->IsDisabled(dg)) continue;
 376       SetCouplingOld(dg,idtrack,h);
 377     }
 378     break;
 379   default:
 380     break;
 381   } // end switch
 382   if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
 383 }
 384
 385 //______________________________________________________________________
 386 void AliITSUSimulationPix::SpreadCharge(Double_t x0,Double_t z0,
 387                                           Int_t ix0,Int_t iz0,
 388                                           Double_t el,Double_t sig,Double_t ld,
 389                                           Int_t t,Int_t hi)
 390 {
 391    // Spreads the charge over neighboring cells. Assume charge is distributed
 392    // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg)
 393    // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig)
 394    // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
 395    // Defined this way, the integral over all x and z is el.
 396    // Inputs:
 397    //    Double_t x0   x position of point where charge is liberated
 398    //    Double_t z0   z position of point where charge is liberated
 399    //    Int_t    ix0  row of cell corresponding to point x0
 400    //    Int_t    iz0  columb of cell corresponding to point z0
 401    //    Double_t el   number of electrons liberated in this step
 402    //    Double_t sig  Sigma difusion for this step (y0 dependent)
 403    //    Double_t ld   lorentz drift in x for this step (y0 dependent)
 404    //    Int_t    t    track number
 405    //    Int_t    ti   hit track index number
 406    //    Int_t    hi   hit "hit" index number
 407    // Outputs:
 408    //     none.
 409    // Return:
 410    //     none.
 411    const Int_t knx = 3,knz = 2;
 412    const Double_t kRoot2 = 1.414213562; // Sqrt(2).
 413    const Double_t kmictocm = 1.0e-4; // convert microns to cm.
 414    Int_t ix,iz,ixs,ixe,izs,ize;
 415    Float_t x,z;  // keep coordinates float (required by AliSegmentation)
 416    Double_t s,sp,x1,x2,z1,z2;
 417    //
 418    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));
 419    if (sig<=0.0) { // if sig<=0 No diffusion to simulate.
 420      UpdateMapSignal(iz0,ix0,t,hi,el);
 421      return;
 422    } // end if
 423    sp = 1.0/(sig*kRoot2);
 424    ixs = TMath::Max(-knx+ix0,0);
 425    ixe = TMath::Min(knx+ix0,fSeg->Npx()-1);
 426    izs = TMath::Max(-knz+iz0,0);
 427    ize = TMath::Min(knz+iz0,fSeg->Npz()-1);
 428    for (ix=ixs;ix<=ixe;ix++)
 429      for (iz=izs;iz<=ize;iz++) {
 430        fSeg->DetToLocal(ix,iz,x,z); // pixel center
 431        double dxi = 0.5*kmictocm*fSeg->Dpx(ix);
 432        double dzi = 0.5*kmictocm*fSeg->Dpz(iz);
 433        x1  = x;
 434        z1  = z;
 435        x2  = x1 + dxi; // Upper
 436        x1 -= dxi;  // Lower
 437        z2  = z1 + dzi; // Upper
 438        z1 -= dzi;  // Lower
 439        x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
 440        x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
 441        z1 -= z0; // Distance from where track traveled
 442        z2 -= z0; // Distance from where track traveled
 443        s   = el*0.25; // Correction based on definision of Erfc
 444        s  *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(sp*x2);
 445        s  *= AliMathBase::ErfcFast(sp*z1) - AliMathBase::ErfcFast(sp*z2);
 446        if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
 447      } // end for ix, iz
 448    //
 449 }
 450
 451 //______________________________________________________________________
 452 void AliITSUSimulationPix::SpreadChargeAsym(Double_t x0,Double_t z0,
 453                                             Int_t ix0,Int_t iz0,
 454                                             Double_t el,Double_t sigx,Double_t sigz,
 455                                             Double_t ld,Int_t t,Int_t hi)
 456 {
 457   // Spreads the charge over neighboring cells. Assume charge is distributed
 458   // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
 459   // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
 460   // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
 461   // Defined this way, the integral over all x and z is el.
 462   // Inputs:
 463   //    Double_t x0   x position of point where charge is liberated
 464   //    Double_t z0   z position of point where charge is liberated
 465   //    Int_t    ix0  row of cell corresponding to point x0
 466   //    Int_t    iz0  columb of cell corresponding to point z0
 467   //    Double_t el   number of electrons liberated in this step
 468   //    Double_t sigx Sigma difusion along x for this step (y0 dependent)
 469   //    Double_t sigz Sigma difusion along z for this step (z0 dependent)
 470   //    Double_t ld   lorentz drift in x for this stip (y0 dependent)
 471   //    Int_t    t    track number
 472   //    Int_t    ti   hit track index number
 473   //    Int_t    hi   hit "hit" index number
 474   // Outputs:
 475   //     none.
 476   // Return:
 477   //     none.
 478   const Int_t knx = 3,knz = 3; // RS: TO TUNE
 479   const Double_t kRoot2 = 1.414213562; // Sqrt(2).
 480   const Double_t kmictocm = 1.0e-4; // convert microns to cm.
 481   Int_t ix,iz,ixs,ixe,izs,ize;
 482   Float_t x,z;   // keep coordinates float (required by AliSegmentation)
 483   Double_t s,spx,spz,x1,x2,z1,z2;
 484   //
 485   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)",
 486                                 x0,z0,ix0,iz0,el,sigx,sigz,t,hi,ld));
 487   if (sigx<=0.0 || sigz<=0.0) { // if sig<=0 No diffusion to simulate.
 488     UpdateMapSignal(iz0,ix0,t,hi,el);
 489     return;
 490   } // end if
 491   spx = 1.0/(sigx*kRoot2);
 492   spz = 1.0/(sigz*kRoot2);
 493   ixs = TMath::Max(-knx+ix0,0);
 494   ixe = TMath::Min(knx+ix0,fSeg->Npx()-1);
 495   izs = TMath::Max(-knz+iz0,0);
 496   ize = TMath::Min(knz+iz0,fSeg->Npz()-1);
 497   for (ix=ixs;ix<=ixe;ix++)
 498     for (iz=izs;iz<=ize;iz++) {
 499       fSeg->DetToLocal(ix,iz,x,z); // pixel center
 500       double dxi = 0.5*kmictocm*fSeg->Dpx(ix);
 501       double dzi = 0.5*kmictocm*fSeg->Dpz(iz);
 502       x1  = x;
 503       z1  = z;
 504       x2  = x1 + dxi; // Upper
 505       x1 -= dxi;      // Lower
 506       z2  = z1 + dzi; // Upper
 507       z1 -= dzi;      // Lower
 508       x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
 509       x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
 510       z1 -= z0; // Distance from where track traveled
 511       z2 -= z0; // Distance from where track traveled
 512       s   = el*0.25; // Correction based on definision of Erfc
 513       s  *= AliMathBase::ErfcFast(spx*x1) - AliMathBase::ErfcFast(spx*x2);
 514       s  *= AliMathBase::ErfcFast(spz*z1) - AliMathBase::ErfcFast(spz*z2);
 515       if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
 516     } // end for ix, iz
 517   //
 518 }
 519
 520 //______________________________________________________________________
 521 void AliITSUSimulationPix::RemoveDeadPixels()
 522 {
 523   // Removes dead pixels on each module (ladder)
 524   // This should be called before going from sdigits to digits (FrompListToDigits)
 525
 526   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibDead();
 527   if (!calObj) return;
 528   //
 529   if (calObj->IsBad()) {ClearMap(); return;} // whole module is masked
 530   //
 531   // remove single bad pixels one by one
 532   int nsingle = calObj->GetNrBadSingle();
 533   UInt_t col,row;
 534   for (int i=nsingle;i--;) {
 535     calObj->GetBadPixelSingle(i,row,col);
 536     fSensMap->DeleteItem(col,row);
 537   }
 538   int nsd = fSensMap->GetEntriesUnsorted();
 539   for (int isd=nsd;isd--;) {
 540     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->AtUnsorted(isd);
 541     if (fSensMap->IsDisabled(sd)) continue;
 542     fSensMap->GetMapIndex(sd->GetUniqueID(),col,row);
 543     int chip = fSeg->GetChipFromChannel(0,col);
 544     //    if (calObj->IsChipMarkedBad(chip)) fSensMap->Disable(sd); // this will simple mark the hit as bad
 545     if (calObj->IsChipMarkedBad(chip)) fSensMap->DeleteItem(sd); // this will suppress hit in the sorted list
 546   }
 547   //
 548 }
 549
 550 //______________________________________________________________________
 551 void AliITSUSimulationPix::AddNoisyPixels()
 552 {
 553   // Adds noisy pixels on each module (ladder)
 554   // This should be called before going from sdigits to digits (FrompListToDigits)
 555   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibNoisy();
 556   if (!calObj) return;
 557   for (Int_t i=calObj->GetNrBad(); i--;) UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i),
 558                                                         10*fSimuParam->GetPixThreshold(fModule));
 559   //
 560 }
 561
 562 //______________________________________________________________________
 563 void AliITSUSimulationPix::FrompListToDigits()
 564 {
 565   // add noise and electronics, perform the zero suppression and add the
 566   // digit to the list
 567   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
 568   UInt_t ix,iz;
 569   Double_t sig;
 570   const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
 571   static AliITSUDigitPix dig;
 572   // RS: in principle:
 573   // 1) for every pixel w/o hit we have to generate a noise and activate the pixel if the noise exceeds the threshold.
 574   // 2) for every pixel with hit we should add random noise and check if the total signal exceeds the threshold
 575   // With many channels this will be too time consuming, hence I do the following
 576   // 1) Precalculate the probability that the nois alone will exceed the threshold.
 577   // 2) Chose randomly empty pixels according to this probability and apply the noise above threshold.
 578   // 3) For pixels having a hits apply the usual noise and compare with threshold
 579   //
 580   // RS may use for ordered random sample generation dl.acm.org/ft_gateway.cfm?id=356313&type=pdf
 581   //
 582   int maxInd = fSensMap->GetMaxIndex();
 583   double minProb = 0.1/maxInd;
 584   //
 585   int nsd = fSensMap->GetEntries();
 586   Int_t prevID=0,curID=0;
 587   TArrayI ordSampleInd(100),ordSample(100);
 588   //
 589   double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(fModule);
 590   fSimuParam->GetPixNoise(fModule, noiseSig, noiseMean);
 591   probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold
 592   //
 593   for (int i=0;i<nsd;i++) {
 594     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
 595     if (fSensMap->IsDisabled(sd)) continue;
 596     curID = (int)sd->GetUniqueID();
 597     //
 598     if (probNoisy>minProb) { // generate randomly noisy pixels above the threshold, with ID's between previous hit and current
 599       CreateNoisyDigits(prevID,curID,probNoisy, noiseSig, noiseMean);
 600       prevID = curID+1;
 601     }
 602     //
 603     if ((sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(fModule)) continue;
 604     if (TMath::Abs(sig)>2147483647.0) { //RS?
 605       //PH 2147483647 is the max. integer
 606       //PH This apparently is a problem which needs investigation
 607       AliWarning(Form("Too big or too small signal value %f",sig));
 608     }
 609     fSensMap->GetMapIndex(sd->GetUniqueID(),iz,ix);
 610     dig.SetCoord1(iz);
 611     dig.SetCoord2(ix);
 612     dig.SetSignal(1);
 613     dig.SetSignalPix((Int_t)sig);
 614     int ntr = sd->GetNTracks();
 615     for (int j=0;j<ntr;j++) {
 616       dig.SetTrack(j,sd->GetTrack(j));
 617       dig.SetHit(j,sd->GetHit(j));
 618     }
 619     for (int j=ntr;j<knmaxtrk;j++) {
 620       dig.SetTrack(j,-3);
 621       dig.SetHit(j,-1);
 622     }
 623     aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix, &dig);
 624   }
 625   // if needed, add noisy pixels with id from last real hit to maxID
 626   if (probNoisy>minProb) CreateNoisyDigits(prevID,maxInd,probNoisy, noiseSig, noiseMean);
 627   //
 628 }
 629
 630 //______________________________________________________________________
 631 Int_t AliITSUSimulationPix::CreateNoisyDigits(Int_t minID,Int_t maxID,double probNoisy, double noise, double base)
 632 {
 633   // create random noisy digits above threshold within id range [minID,maxID[
 634   // see FrompListToDigits for details
 635   //
 636   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
 637   UInt_t ix,iz;
 638   static AliITSUDigitPix dig;
 639   static TArrayI ordSampleInd(100),ordSample(100); //RS!!! static is not thread-safe!!!
 640   const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
 641   //
 642   Int_t ncand = 0;
 643   int npix = maxID-minID;
 644   if (npix<1 || (ncand=gRandom->Poisson(npix*probNoisy))<1) return ncand; // decide how many noisy pixels will be added
 645   ncand = GenOrderedSample(npix,ncand,ordSample,ordSampleInd);
 646   int* ordV = ordSample.GetArray();
 647   int* ordI = ordSampleInd.GetArray();
 648   for (int j=0;j<ncand;j++) {
 649     fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],iz,ix);   // create noisy digit
 650     dig.SetCoord1(iz);
 651     dig.SetCoord2(ix);
 652     dig.SetSignal(1);
 653     dig.SetSignalPix((Int_t)AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,base,noise));
 654     for (int k=knmaxtrk;k--;) {
 655       dig.SetTrack(k,-3);
 656       dig.SetHit(k,-1);
 657     }
 658     aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix,&dig);
 659     if (GetDebug(2)) AliInfo(Form("Add noisy pixel %d(%d/%d) Noise=%d",ordV[ordI[j]],iz,ix,dig.GetSignalPix()));
 660   }
 661   return ncand;
 662 }
 663
 664 //______________________________________________________________________
 665 void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old, Int_t ntrack, Int_t idhit)
 666 {
 667   //  Take into account the coupling between adiacent pixels.
 668   //  The parameters probcol and probrow are the probability of the
 669   //  signal in one pixel shared in the two adjacent pixels along
 670   //  the column and row direction, respectively.
 671   //  Note pList is goten via GetMap() and module is not need any more.
 672   //  Otherwise it is identical to that coded by Tiziano Virgili (BSN).
 673   //Begin_Html
 674   /*
 675     <img src="picts/ITS/barimodel_3.gif">
 676      </pre>
 677      <br clear=left>
 678      <font size=+2 color=red>
 679      <a href="mailto:tiziano.virgili@cern.ch"></a>.
 680      </font>
 681      <pre>
 682    */
 683    //End_Html
 684    // Inputs:
 685   // old                  existing AliITSUSDigit
 686   // Int_t ntrack         track incex number
 687   // Int_t idhit          hit index number
 688   UInt_t col,row;
 689   Double_t pulse1,pulse2;
 690   Double_t couplR=0.0,couplC=0.0;
 691   Double_t xr=0.;
 692   //
 693   fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
 694   fSimuParam->GetPixCouplingParam(couplC,couplR);
 695   if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
 696                                 col,row,ntrack,idhit,couplC,couplR));
 697   pulse2 = pulse1 = old->GetSignal();
 698   if (pulse1<fSimuParam->GetPixMinElToAdd()) return; // too small signal
 699   for (Int_t isign=-1;isign<=1;isign+=2) {
 700     //
 701     // loop in col direction
 702     int j1 = int(col) + isign;
 703     xr = gRandom->Rndm();
 704     if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
 705     //
 706     // loop in row direction
 707     int j2 = int(row) + isign;
 708     xr = gRandom->Rndm();
 709     if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
 710   }
 711   //
 712 }
 713
 714 //______________________________________________________________________
 715 void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, Int_t ntrack,Int_t idhit)
 716 {
 717   //  Take into account the coupling between adiacent pixels.
 718   //  The parameters probcol and probrow are the fractions of the
 719   //  signal in one pixel shared in the two adjacent pixels along
 720   //  the column and row direction, respectively.
 721   //Begin_Html
 722   /*
 723     <img src="picts/ITS/barimodel_3.gif">
 724     </pre>
 725     <br clear=left>
 726     <font size=+2 color=red>
 727     <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
 728     </font>
 729     <pre>
 730   */
 731   //End_Html
 732   // Inputs:
 733   // old            existing AliITSUSDigit
 734   // ntrack         track incex number
 735   // idhit          hit index number
 736   // module         module number
 737   //
 738   UInt_t col,row;
 739   Double_t pulse1,pulse2;
 740   Double_t couplR=0.0,couplC=0.0;
 741   //
 742   fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
 743   fSimuParam->GetPixCouplingParam(couplC,couplR);
 744   if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
 745                                 col,row,ntrack,idhit,couplC,couplR));
 746  //
 747  if (old->GetSignal()<fSimuParam->GetPixMinElToAdd()) return; // too small signal
 748  for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
 749    pulse2 = pulse1 = old->GetSignal();
 750    //
 751    int j1 = int(col)+isign;
 752    pulse1 *= couplC;
 753    if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(fModule))) pulse1 = old->GetSignal();
 754    else UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
 755
 756    // loop in row direction
 757    int j2 = int(row) + isign;
 758    pulse2 *= couplR;
 759    if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(fModule))) pulse2 = old->GetSignal();
 760    else UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
 761  } // for isign
 762 }
 763
 764 //______________________________________________________________________
 765 void AliITSUSimulationPix::GenerateStrobePhase()
 766 {
 767   // Generate randomly the strobe
 768   // phase w.r.t to the LHC clock
 769   // Done once per event
 770   const Double_t kBunchLenght = 25e-9; // LHC clock
 771   fStrobePhase = ((Double_t)gRandom->Integer(fStrobeLenght))*kBunchLenght-
 772     (Double_t)fStrobeLenght*kBunchLenght+
 773     kBunchLenght/2;
 774 }
 775