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