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 <TH2.h>
  25 #include <TString.h>
  26 #include "AliITSU.h"
  27 #include "AliITSUDigitPix.h"
  28 #include "AliITSUHit.h"
  29 #include "AliITSUModule.h"
  30 #include "AliITSUSensMap.h"
  31 #include "AliITSUCalibrationPix.h"
  32 #include "AliITSUSegmentationPix.h"
  33 #include "AliITSUSimulationPix.h"
  34 #include "AliLog.h"
  35 #include "AliRun.h"
  36 #include "AliMagF.h"
  37 #include "AliMathBase.h"
  38 #include "AliITSUSimuParam.h"
  39 #include "AliITSUSDigit.h"
  40 #include "AliITSUParamList.h"
  41
  42 using std::cout;
  43 using std::endl;
  44 using namespace TMath;
  45
  46 ClassImp(AliITSUSimulationPix)
  47 ////////////////////////////////////////////////////////////////////////
  48 //  Version: 1
  49 //  Modified by D. Elia, G.E. Bruno, H. Tydesjo
  50 //  Fast diffusion code by Bjorn S. Nilsen
  51 //  March-April 2006
  52 //  October     2007: GetCalibrationObjects() removed
  53 //
  54 //  Version: 0
  55 //  Written by Boris Batyunya
  56 //  December 20 1999
  57 //
  58 //  Adapted for pixels of ITS upgrade July 2012, ruben.shahoyan@cern.ch
  59 //
  60 //  AliITSUSimulationPix is to do the simulation of pixels
  61 //
  62 //  2013 Feb: Added MonoPix response and nois calculation al la MIMOSA32 (levente.molnar@cern.ch)
  63 //
  64 ////////////////////////////////////////////////////////////////////////
  65
  66 //______________________________________________________________________
  67 AliITSUSimulationPix::AliITSUSimulationPix()
  68 :  fTanLorAng(0)
  69   ,fReadOutCycleLength(25e-6)
  70   ,fReadOutCycleOffset(0)
  71   ,fGlobalChargeScale(1.0)
  72   ,fSpread2DHisto(0)
  73   ,fSpreadFun(0)
  74   ,fROTimeFun(0)
  75 {
  76    // Default constructor.
  77   SetUniqueID(AliITSUGeomTGeo::kDetTypePix);
  78 }
  79
  80 //______________________________________________________________________
  81 AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap* map)
  82   :AliITSUSimulation(sim,map)
  83   ,fTanLorAng(0)
  84   ,fReadOutCycleLength(25e-6)
  85   ,fReadOutCycleOffset(0)
  86   ,fGlobalChargeScale(1.0)
  87   ,fSpread2DHisto(0)
  88   ,fSpreadFun(0)
  89   ,fROTimeFun(0)
  90 {
  91   // standard constructor
  92   SetUniqueID(AliITSUGeomTGeo::kDetTypePix);
  93   Init();
  94 }
  95
  96 //______________________________________________________________________
  97 AliITSUSimulationPix::AliITSUSimulationPix(const AliITSUSimulationPix &s)
  98   :AliITSUSimulation(s)
  99   ,fTanLorAng(s.fTanLorAng)
 100   ,fReadOutCycleLength(s.fReadOutCycleLength)
 101   ,fReadOutCycleOffset(s.fReadOutCycleOffset)
 102   ,fGlobalChargeScale(s.fGlobalChargeScale)
 103   ,fSpread2DHisto(s.fSpread2DHisto)
 104   ,fSpreadFun(s.fSpreadFun)
 105   ,fROTimeFun(s.fROTimeFun)
 106 {
 107   //     Copy Constructor
 108 }
 109
 110
 111 //______________________________________________________________________
 112 AliITSUSimulationPix::~AliITSUSimulationPix()
 113 {
 114   // destructor
 115   // only the sens map is owned and it is deleted by ~AliITSUSimulation
 116 }
 117
 118 //______________________________________________________________________
 119 AliITSUSimulationPix& AliITSUSimulationPix::operator=(const AliITSUSimulationPix &s)
 120 {
 121   //    Assignment operator
 122   if (&s == this) return *this;
 123   AliITSUSimulation::operator=(s);
 124   fReadOutCycleLength = s.fReadOutCycleLength;
 125   fReadOutCycleOffset = s.fReadOutCycleOffset;
 126   fSpread2DHisto = s.fSpread2DHisto;
 127   fGlobalChargeScale = s.fGlobalChargeScale;
 128   fSpreadFun    = s.fSpreadFun;
 129   fROTimeFun    = s.fROTimeFun;
 130   //
 131   return *this;
 132 }
 133
 134 //______________________________________________________________________
 135 void AliITSUSimulationPix::Init()
 136 {
 137   // Initilization
 138   if (fSimuParam->GetPixLorentzDrift()) SetTanLorAngle(fSimuParam->GetPixLorentzHoleWeight());
 139 }
 140
 141 //______________________________________________________________________
 142 Bool_t AliITSUSimulationPix::SetTanLorAngle(Double_t weightHole)
 143 {
 144   // This function set the Tangent of the Lorentz angle.
 145   // A weighted average is used for electrons and holes
 146   // Input: Double_t weightHole: wheight for hole: it should be in the range [0,1]
 147   // output: Bool_t : kTRUE in case of success
 148   //
 149   if (weightHole<0) {
 150     weightHole=0.;
 151     AliWarning("You have asked for negative Hole weight");
 152     AliWarning("I'm going to use only electrons");
 153   }
 154   if (weightHole>1) {
 155     weightHole=1.;
 156     AliWarning("You have asked for weight > 1");
 157     AliWarning("I'm going to use only holes");
 158   }
 159   Double_t weightEle=1.-weightHole;
 160   AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 161   if (!fld) AliFatal("The field is not initialized");
 162   Double_t bz = fld->SolenoidField();
 163   fTanLorAng = Tan(weightHole*fSimuParam->LorentzAngleHole(bz) +
 164                           weightEle*fSimuParam->LorentzAngleElectron(bz));
 165   return kTRUE;
 166 }
 167
 168 //_____________________________________________________________________
 169 void AliITSUSimulationPix::SDigitiseModule()
 170 {
 171   //  This function begins the work of creating S-Digits.
 172
 173   AliDebug(10,Form("In event %d module %d there are %d hits", fEvent, fModule->GetIndex(),fModule->GetNHits()));
 174   //
 175   if (fModule->GetNHits()) Hits2SDigitsFast();
 176   //
 177   if (fSimuParam->GetPixAddNoisyFlag())   AddNoisyPixels();
 178   if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();
 179   if (!fSensMap->GetEntries()) return;
 180   WriteSDigits();
 181   ClearMap();
 182   //
 183 }
 184
 185 //______________________________________________________________________
 186 void AliITSUSimulationPix::WriteSDigits()
 187 {
 188   //  This function adds each S-Digit to pList
 189   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
 190   int nsd = fSensMap->GetEntries();
 191   for (int i=0;i<nsd;i++) {
 192     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
 193     if (!sd->GetSumSignal()>0 || fSensMap->IsDisabled(sd)) continue;
 194     aliITS->AddSumDigit(*sd);
 195   }
 196   return;
 197 }
 198
 199 //______________________________________________________________________
 200 void AliITSUSimulationPix::FinishSDigitiseModule()
 201 {
 202    //  This function calls SDigitsToDigits which creates Digits from SDigits
 203   FrompListToDigits();
 204   ClearMap();
 205   return;
 206 }
 207
 208 //______________________________________________________________________
 209 void AliITSUSimulationPix::DigitiseModule()
 210 {
 211   //  This function creates Digits straight from the hits and then adds
 212   //  electronic noise to the digits before adding them to pList
 213   //  Each of the input variables is passed along to Hits2SDigits
 214   //
 215   // pick charge spread function
 216   Hits2SDigitsFast();
 217   //
 218   if (fSimuParam->GetPixAddNoisyFlag())   AddNoisyPixels();
 219   if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();
 220   FrompListToDigits();
 221   ClearMap();
 222 }
 223
 224 //______________________________________________________________________
 225 void AliITSUSimulationPix::Hits2SDigits()
 226 {
 227   // Does the charge distributions using Gaussian diffusion charge charing.
 228   Int_t nhits = fModule->GetNHits();
 229   if (!nhits) return;
 230   //
 231   Int_t h,ix,iz,i;
 232   Int_t idtrack;
 233   Float_t x,y,z;  // keep coordinates float (required by AliSegmentation)
 234   Double_t tof,x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0;
 235   Double_t t,tp,st,dt=0.2,el;
 236   Double_t thick = 0.5*fSeg->Dy();  // Half Thickness
 237
 238   //
 239   for (h=0;h<nhits;h++) {
 240     //
 241     if (!fModule->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
 242     st = Sqrt(x1*x1+y1*y1+z1*z1);
 243     if (st>0.0) {
 244       st = (Double_t)((Int_t)(st*1e4)); // number of microns
 245       if (st<=1.0) st = 1.0;
 246       dt = 1.0/st;               // RS TODO: do we need 1 micron steps?
 247       double dy = dt*thick;
 248       y = -0.5*dy;
 249       for (t=0.0;t<1.0;t+=dt) { // Integrate over t
 250         tp  = t+0.5*dt;
 251         x   = x0+x1*tp;
 252         y  += dy;
 253         z   = z0+z1*tp;
 254         if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 255         el  = dt * de / fSimuParam->GetGeVToCharge();
 256         //
 257         if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
 258         SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
 259       } // end for t
 260     } else { // st == 0.0 deposit it at this point
 261       x   = x0;
 262       y   = y0 + 0.5*thick;
 263       z   = z0;
 264       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 265       el  = de / fSimuParam->GetGeVToCharge();
 266       if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
 267       SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
 268     } // end if st>0.0
 269   } // Loop over all hits h
 270   //
 271   // Coupling
 272   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
 273   AliITSUSDigit* dg = 0;
 274   switch (fSimuParam->GetPixCouplingOption()) {
 275   case AliITSUSimuParam::kNewCouplingPix :
 276     for (i=nd;i--;) {
 277       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 278       if (fSensMap->IsDisabled(dg)) continue;
 279       SetCoupling(dg,idtrack,h);
 280     }
 281     break;
 282   case AliITSUSimuParam::kOldCouplingPix:
 283     for (i=nd;i--;) {
 284       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 285       if (fSensMap->IsDisabled(dg)) continue;
 286       SetCouplingOld(dg,idtrack,h);
 287     }
 288     break;
 289   default:
 290     break;
 291
 292   } // end switch
 293   if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
 294 }
 295
 296 //______________________________________________________________________
 297 void AliITSUSimulationPix::Hits2SDigitsFast()
 298 {
 299   // Does the charge distributions using Gaussian diffusion charge charing.    // Inputs:
 300   //    AliITSUModule *mod  Pointer to this module
 301   //
 302   TObjArray *hits = fModule->GetHits();
 303   Int_t nhits = hits->GetEntriesFast();
 304   if (nhits<=0) return;
 305   //
 306   Int_t h,ix,iz,i;
 307   Int_t idtrack;
 308   Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
 309   Double_t tof,x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0;
 310   Double_t step,st,el,de=0.0;
 311   Double_t minDim = Min(fSeg->Dpx(1),fSeg->Dpz(1)); // RStmp: smallest pitch
 312   Double_t thick = fSeg->Dy();
 313   //
 314   for (h=0;h<nhits;h++) {
 315     //
 316     if (!fModule->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
 317     //
 318     st = Sqrt(x1*x1+y1*y1+z1*z1);
 319     if (st>0.0) {
 320       int np = int(1.5*st/minDim);  //RStmp: inject the points in such a way that there is ~1.5 point per cell
 321       np = TMath::Max(1.0*np,fResponseParam->GetParameter(kSpreadFunMinSteps));
 322       AliDebug(10,Form(" Number of charge injection steps is set to %d ",np));
 323       double dstep = 1./np;
 324       double dy = dstep*thick;
 325       y = -0.5*dy;
 326       step = -0.5*dstep;
 327       for (i=0;i<np;i++) {          //RStmp Integrate over t
 328         //      for (i=0;i<kn10;i++) { // Integrate over t
 329         step  += dstep;  // RStmp kti[i];
 330         x   = x0+x1*step;
 331         y  += dy;
 332         z   = z0+z1*step;
 333         if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 334         el  = fGlobalChargeScale * dstep * de/fSimuParam->GetGeVToCharge();
 335         if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
 336         SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
 337       } // end for i // End Integrate over t
 338     }
 339     else { // st == 0.0 deposit it at this point
 340       x   = x0;
 341       y   = y0+0.5*thick;
 342       z   = z0;
 343       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 344       el  = de / fSimuParam->GetGeVToCharge();
 345       if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
 346       SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
 347     } // end if st>0.0
 348
 349   } // Loop over all hits h
 350
 351   // Coupling
 352   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
 353   AliITSUSDigit* dg = 0;
 354   switch (fSimuParam->GetPixCouplingOption()) {
 355   case AliITSUSimuParam::kNewCouplingPix :
 356     for (i=nd;i--;) {
 357       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 358       if (fSensMap->IsDisabled(dg)) continue;
 359       SetCoupling(dg,idtrack,h);
 360     }
 361   case AliITSUSimuParam::kOldCouplingPix:
 362     for (i=nd;i--;) {
 363       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 364       if (fSensMap->IsDisabled(dg)) continue;
 365       SetCouplingOld(dg,idtrack,h);
 366     }
 367     break;
 368   default:
 369     break;
 370   } // end switch
 371   if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
 372 }
 373
 374 //______________________________________________________________________
 375 void AliITSUSimulationPix::SpreadCharge2D(Double_t x0,Double_t z0, Double_t dy, Int_t ix0,Int_t iz0,
 376                                           Double_t el, Double_t tof, Int_t tID, Int_t hID)
 377 {
 378   // Spreads the charge over neighboring cells. Assume charge is distributed
 379   // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
 380   // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
 381   // Defined this way, the integral over all x and z is el.
 382   // Inputs:
 383   //    Double_t x0   x position of point where charge is liberated (local)
 384   //    Double_t z0   z position of point where charge is liberated (local)
 385   //    Double_t dy   distance from the entrance surface (diffusion sigma may depend on it)
 386   //    Int_t    ix0  row of cell corresponding to point x0
 387   //    Int_t    iz0  columb of cell corresponding to point z0
 388   //    Double_t el   number of electrons liberated in this step
 389   //    Double_t sigx Sigma difusion along x for this step (y0 dependent)
 390   //    Double_t sigz Sigma difusion along z for this step (z0 dependent)
 391   //    Int_t    tID  track number
 392   //    Int_t    hID  hit "hit" index number
 393   //
 394   Int_t ix,iz,ixs,ixe,izs,ize;
 395   Float_t x,z;   // keep coordinates float (required by AliSegmentation)
 396   Float_t xdioshift = 0 , zdioshift = 0 ;
 397   Double_t s,dtIn[kNDtSpread]; // data transfered to spread function for integral calculation
 398   //
 399   if (GetDebug(2)) AliInfo(Form("(x0=%e,z0=%e,dy=%e, ix0=%d,iz0=%d,el=%e,tID=%d,hID=%d)",x0,z0,dy,ix0,iz0,el,tID,hID));
 400   //
 401   Double_t &x1 = dtIn[kCellX1];
 402   Double_t &x2 = dtIn[kCellX2];
 403   Double_t &z1 = dtIn[kCellZ1];
 404   Double_t &z2 = dtIn[kCellZ2];
 405   //
 406   int nx = GetResponseParam()->GetParameter(kSpreadFunParamNXoffs);
 407   int nz = GetResponseParam()->GetParameter(kSpreadFunParamNZoffs);
 408   //
 409   dtIn[kCellYDepth]  = dy;
 410   ixs = Max(-nx+ix0,0);
 411   ixe = Min( nx+ix0,fSeg->Npx()-1);
 412   izs = Max(-nz+iz0,0);
 413   ize = Min( nz+iz0,fSeg->Npz()-1);
 414   for (ix=ixs;ix<=ixe;ix++)
 415     for (iz=izs;iz<=ize;iz++) {
 416       //
 417       // Check if the hit is inside readout window
 418       int cycleRO = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof);
 419       if (Abs(cycleRO)>kMaxROCycleAccept) continue;
 420       //
 421       fSeg->DetToLocal(ix,iz,x,z); // pixel center
 422       xdioshift = zdioshift = 0;
 423       CalcDiodeShiftInPixel(ix,iz,xdioshift,zdioshift);    // Check and apply diode shift if needed
 424       double dxi = 0.5*fSeg->Dpx(ix);
 425       double dzi = 0.5*fSeg->Dpz(iz);
 426       x1  = (x + xdioshift) - x0;   // calculate distance of cell boundaries from injection center
 427       z1  = (z + zdioshift) - z0;
 428       x2  = x1 + dxi; // Upper
 429       x1 -= dxi;      // Lower
 430       z2  = z1 + dzi; // Upper
 431       z1 -= dzi;      // Lower
 432       s   = el* (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fSpreadFun)(dtIn);
 433       if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,tID,hID,s,cycleRO);
 434     } // end for ix, iz
 435   //
 436 }
 437
 438 //______________________________________________________________________
 439 Double_t AliITSUSimulationPix::SpreadFunDoubleGauss2D(const Double_t *dtIn)
 440 {
 441   // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
 442   // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
 443   // The spread function is assumed to be double gaussian in 2D
 444   // Parameters should be: mean0,sigma0, mean1,sigma1, relative strenght of 2nd gaussian wrt 1st one
 445   //
 446   // 1st gaussian
 447   double intg1 = GausInt2D(fResponseParam->GetParameter(kG2SigX0),  // sigX
 448                            dtIn[kCellX1]-fResponseParam->GetParameter(kG2MeanX0),      // x1-xmean
 449                            dtIn[kCellX2]-fResponseParam->GetParameter(kG2MeanX0),      // x2-xmean
 450                            fResponseParam->GetParameter(kG2SigZ0),  // sigZ
 451                            dtIn[kCellZ1]-fResponseParam->GetParameter(kG2MeanZ0),    // z1-zmean
 452                            dtIn[kCellZ2]-fResponseParam->GetParameter(kG2MeanZ0));   // z2-zmean
 453   // 2nd gaussian
 454   double intg2 = GausInt2D(fResponseParam->GetParameter(kG2SigX1),  // sigX
 455                            dtIn[kCellX1]-fResponseParam->GetParameter(kG2MeanX1),      // x1-xmean
 456                            dtIn[kCellX2]-fResponseParam->GetParameter(kG2MeanX1),      // x2-xmean
 457                            fResponseParam->GetParameter(kG2SigZ1),  // sigZ
 458                            dtIn[kCellZ1]-fResponseParam->GetParameter(kG2MeanZ1),    // z1-zmean
 459                            dtIn[kCellZ2]-fResponseParam->GetParameter(kG2MeanZ1));   // z2-zmean
 460   double scl = fResponseParam->GetParameter(kG2ScaleG2);
 461   return (intg1+intg2*scl)/(1+scl);
 462   //
 463 }
 464
 465
 466 //______________________________________________________________________
 467 Double_t AliITSUSimulationPix::SpreadFrom2DHisto(const Double_t *dtIn)
 468 {
 469   // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
 470   // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
 471   // The spread function integral is taken from fSpread2DHisto extracted from the sensor response parameters
 472   // list in the method SetResponseParam. The histo must return the fraction of charge integrates in the
 473   // cell whose center is on the distance X=(dtIn[kCellX1]+dtIn[kCellX2])/2 and Z=(dtIn[kCellZ1]+dtIn[kCellZ2])/2
 474   // from the injection point.
 475   //
 476   Double_t qpixfrac = 0;
 477   Double_t xintp = 1e4*(dtIn[kCellX1]+dtIn[kCellX2])/2.0;
 478   Double_t zintp = 1e4*(dtIn[kCellZ1]+dtIn[kCellZ2])/2.0;
 479   //
 480   qpixfrac =  fSpread2DHisto->Interpolate(xintp,zintp); //the PSF map is given in um but the dtIn is in cm so we need to convert it
 481   //
 482   return qpixfrac;
 483 }
 484
 485 //______________________________________________________________________
 486 Double_t AliITSUSimulationPix::SpreadFunGauss2D(const Double_t *dtIn)
 487 {
 488   // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
 489   // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
 490   // The spread function is assumed to be gaussian in 2D
 491   // Parameters should be: mean0,sigma0
 492   return GausInt2D(fResponseParam->GetParameter(kG1SigX),  // sigX
 493                    fResponseParam->GetParameter(kG1SigZ),  // sigZ
 494                    dtIn[kCellX1]-fResponseParam->GetParameter(kG1MeanX),    // x1-xmean
 495                    dtIn[kCellX2]-fResponseParam->GetParameter(kG1MeanX),    // x2-xmean
 496                    dtIn[kCellZ1]-fResponseParam->GetParameter(kG1MeanZ),    // z1-zmean
 497                    dtIn[kCellZ2]-fResponseParam->GetParameter(kG1MeanZ));   // z2-zmean
 498   //
 499 }
 500
 501 //______________________________________________________________________
 502 void AliITSUSimulationPix::RemoveDeadPixels()
 503 {
 504   // Removes dead pixels on each module (ladder)
 505   // This should be called before going from sdigits to digits (FrompListToDigits)
 506
 507   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibDead();
 508   if (!calObj) return;
 509   //
 510   if (calObj->IsBad()) {ClearMap(); return;} // whole module is masked
 511   //
 512   // remove single bad pixels one by one
 513   int nsingle = calObj->GetNrBadSingle();
 514   UInt_t col,row,cycle;
 515   for (int i=nsingle;i--;) {
 516     calObj->GetBadPixelSingle(i,row,col);
 517     fSensMap->DeleteItem(col,row);
 518   }
 519   int nsd = fSensMap->GetEntriesUnsorted();
 520   for (int isd=nsd;isd--;) {
 521     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->AtUnsorted(isd);
 522     if (fSensMap->IsDisabled(sd)) continue;
 523     fSensMap->GetMapIndex(sd->GetUniqueID(),col,row,cycle);
 524     int chip = fSeg->GetChipFromChannel(0,col);
 525     //    if (calObj->IsChipMarkedBad(chip)) fSensMap->Disable(sd); // this will simple mark the hit as bad
 526     if (calObj->IsChipMarkedBad(chip)) fSensMap->DeleteItem(sd); // this will suppress hit in the sorted list
 527   }
 528   //
 529 }
 530
 531 //______________________________________________________________________
 532 void AliITSUSimulationPix::AddNoisyPixels()
 533 {
 534   // Adds noisy pixels on each module (ladder)
 535   // This should be called before going from sdigits to digits (FrompListToDigits)
 536   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibNoisy();
 537   if (!calObj) { AliDebug(10,Form("  No Calib Object for Noise!!! ")); return;}
 538   for (Int_t i=calObj->GetNrBad(); i--;) UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i),
 539                                                         10*fSimuParam->GetPixThreshold(fModule->GetIndex()));
 540   //
 541 }
 542
 543 //______________________________________________________________________
 544 void AliITSUSimulationPix::FrompListToDigits()
 545 {
 546   // add noise and electronics, perform the zero suppression and add the
 547   // digit to the list
 548   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
 549   UInt_t ix,iz,iCycle;
 550   Double_t sig;
 551   const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
 552   static AliITSUDigitPix dig;
 553   // RS: in principle:
 554   // 1) for every pixel w/o hit we have to generate a noise and activate the pixel if the noise exceeds the threshold.
 555   // 2) for every pixel with hit we should add random noise and check if the total signal exceeds the threshold
 556   // With many channels this will be too time consuming, hence I do the following
 557   // 1) Precalculate the probability that the nois alone will exceed the threshold (probnoisy).
 558   // 2) Chose randomly empty pixels according to fakes rate
 559   // 3) For pixels having a hits apply the usual noise and compare with threshold
 560   //
 561   // RS may use for ordered random sample generation dl.acm.org/ft_gateway.cfm?id=356313&type=pdf
 562   //
 563   int maxInd = fSensMap->GetMaxIndex();
 564   int modId = fModule->GetIndex();
 565   //
 566   int nsd = fSensMap->GetEntries();
 567   Int_t prevID=0,curID=0;
 568   TArrayI ordSampleInd(100),ordSample(100);
 569   //
 570   double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(modId);
 571   fSimuParam->GetPixNoise(modId, noiseSig, noiseMean);
 572   probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold
 573   //
 574   for (int i=0;i<nsd;i++) {
 575     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
 576     if (fSensMap->IsDisabled(sd)) continue;
 577     curID = (int)sd->GetUniqueID();
 578     //
 579     if ( fSimuParam->GetPixAddNoisyFlag() ) {
 580       CreateNoisyDigits(prevID,curID,probNoisy, noiseSig, noiseMean);
 581       prevID = curID+1;
 582       //
 583       // add noise also to sdigit with signal
 584       sd->AddNoise(AliITSUSimuParam::GenerateNoiseQFunction(0,noiseMean,noiseSig));
 585     }
 586     //
 587     if ((sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(modId)) continue;
 588     if (Abs(sig)>2147483647.0) { //RS?
 589       //PH 2147483647 is the max. integer
 590       //PH This apparently is a problem which needs investigation
 591       AliWarning(Form("Too big or too small signal value %f",sig));
 592     }
 593     fSensMap->GetMapIndex(sd->GetUniqueID(),iz,ix,iCycle);
 594     dig.SetCoord1(iz);
 595     dig.SetCoord2(ix);
 596     dig.SetROCycle(int(iCycle) - kMaxROCycleAccept);
 597     dig.SetSignal((Int_t)sig);
 598     dig.SetSignalPix((Int_t)sig);
 599     int ntr = sd->GetNTracks();
 600     for (int j=0;j<ntr;j++) {
 601       dig.SetTrack(j,sd->GetTrack(j));
 602       dig.SetHit(j,sd->GetHit(j));
 603     }
 604     for (int j=ntr;j<knmaxtrk;j++) {
 605       dig.SetTrack(j,-3);
 606       dig.SetHit(j,-1);
 607     }
 608     aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix, &dig);
 609   }
 610   // if needed, add noisy pixels with id from last real hit to maxID
 611   if (fSimuParam->GetPixAddNoisyFlag() &&
 612       (nsd>0  || ( nsd==0 && fSimuParam->GetPixNoiseInAllMod()) )) {
 613     AliDebug(10,Form("CreateNoisyDigits2( %d ,%d) module %d",prevID,maxInd,modId));
 614     CreateNoisyDigits(prevID,maxInd,probNoisy, noiseSig, noiseMean);
 615   }
 616   //
 617 }
 618
 619 //______________________________________________________________________
 620 Int_t AliITSUSimulationPix::CreateNoisyDigits(Int_t minID,Int_t maxID,double probNoisy, double noise, double base)
 621 {
 622   // create random noisy digits above threshold within id range [minID,maxID[
 623   // see FrompListToDigits for details
 624   //
 625   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
 626   UInt_t ix,iz,iCycle;
 627   static AliITSUDigitPix dig;
 628   static TArrayI ordSampleInd(100),ordSample(100); //RS!!! static is not thread-safe!!!
 629   const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
 630   //
 631   Int_t ncand = 0;
 632   int npix = maxID-minID;
 633   if (npix<1 || (ncand=gRandom->Poisson(npix * fSimuParam->GetPixFakeRate() ))<1) return ncand; // decide how many noisy pixels will be added
 634   ncand = GenOrderedSample(npix,ncand,ordSample,ordSampleInd);
 635   int* ordV = ordSample.GetArray();
 636   int* ordI = ordSampleInd.GetArray();
 637   for (int j=0;j<ncand;j++) {
 638     fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],iz,ix,iCycle);   // create noisy digit
 639     dig.SetCoord1(iz);
 640     dig.SetCoord2(ix);
 641     dig.SetROCycle(int(iCycle) - kMaxROCycleAccept);
 642     dig.SetSignal(1);
 643     dig.SetSignalPix((Int_t)AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,base,noise));
 644     for (int k=knmaxtrk;k--;) {
 645       dig.SetTrack(k,-3);
 646       dig.SetHit(k,-1);
 647     }
 648     aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix,&dig);
 649     AliDebug(10,Form("Add noisy pixel %d(%d/%d) Noise=%d",ordV[ordI[j]],iz,ix,dig.GetSignalPix()));
 650   }
 651   return ncand;
 652 }
 653
 654 //______________________________________________________________________
 655 void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old, Int_t ntrack, Int_t idhit)
 656 {
 657   //  Take into account the coupling between adiacent pixels.
 658   //  The parameters probcol and probrow are the probability of the
 659   //  signal in one pixel shared in the two adjacent pixels along
 660   //  the column and row direction, respectively.
 661   //  Note pList is goten via GetMap() and module is not need any more.
 662   //  Otherwise it is identical to that coded by Tiziano Virgili (BSN).
 663   //Begin_Html
 664   /*
 665     <img src="picts/ITS/barimodel_3.gif">
 666      </pre>
 667      <br clear=left>
 668      <font size=+2 color=red>
 669      <a href="mailto:tiziano.virgili@cern.ch"></a>.
 670      </font>
 671      <pre>
 672    */
 673    //End_Html
 674    // Inputs:
 675   // old                  existing AliITSUSDigit
 676   // Int_t ntrack         track incex number
 677   // Int_t idhit          hit index number
 678   UInt_t col,row,iCycle;
 679   Double_t pulse1,pulse2;
 680   Double_t couplR=0.0,couplC=0.0;
 681   Double_t xr=0.;
 682   //
 683   fSensMap->GetMapIndex(old->GetUniqueID(),col,row,iCycle);
 684   int cycle = int(iCycle)-kMaxROCycleAccept;
 685   fSimuParam->GetPixCouplingParam(couplC,couplR);
 686   if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
 687                                 col,row,ntrack,idhit,couplC,couplR));
 688   pulse2 = pulse1 = old->GetSignal();
 689   if (pulse1<fSimuParam->GetPixMinElToAdd()) return; // too small signal
 690   for (Int_t isign=-1;isign<=1;isign+=2) {
 691     //
 692     // loop in col direction
 693     int j1 = int(col) + isign;
 694     xr = gRandom->Rndm();
 695     if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1,cycle);
 696     //
 697     // loop in row direction
 698     int j2 = int(row) + isign;
 699     xr = gRandom->Rndm();
 700     if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2,cycle);
 701   }
 702   //
 703 }
 704
 705 //______________________________________________________________________
 706 void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, Int_t ntrack,Int_t idhit)
 707 {
 708   //  Take into account the coupling between adiacent pixels.
 709   //  The parameters probcol and probrow are the fractions of the
 710   //  signal in one pixel shared in the two adjacent pixels along
 711   //  the column and row direction, respectively.
 712   //Begin_Html
 713   /*
 714     <img src="picts/ITS/barimodel_3.gif">
 715     </pre>
 716     <br clear=left>
 717     <font size=+2 color=red>
 718     <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
 719     </font>
 720     <pre>
 721   */
 722   //End_Html
 723   // Inputs:
 724   // old            existing AliITSUSDigit
 725   // ntrack         track incex number
 726   // idhit          hit index number
 727   // module         module number
 728   //
 729   UInt_t col,row,iCycle;
 730   Int_t modId = fModule->GetIndex();
 731   Double_t pulse1,pulse2;
 732   Double_t couplR=0.0,couplC=0.0;
 733   //
 734   fSensMap->GetMapIndex(old->GetUniqueID(),col,row,iCycle);
 735   int cycle = int(iCycle)-kMaxROCycleAccept;
 736   fSimuParam->GetPixCouplingParam(couplC,couplR);
 737   if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,roCycle=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
 738                                 col,row,iCycle,ntrack,idhit,couplC,couplR));
 739  //
 740  if (old->GetSignal()<fSimuParam->GetPixMinElToAdd()) return; // too small signal
 741  for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
 742    pulse2 = pulse1 = old->GetSignal();
 743    //
 744    int j1 = int(col)+isign;
 745    pulse1 *= couplC;
 746    if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(modId))) pulse1 = old->GetSignal();
 747    else UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1,cycle);
 748
 749    // loop in row direction
 750    int j2 = int(row) + isign;
 751    pulse2 *= couplR;
 752    if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(modId))) pulse2 = old->GetSignal();
 753    else UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2,cycle);
 754  } // for isign
 755 }
 756
 757 //______________________________________________________________________
 758 void AliITSUSimulationPix::GenerateReadOutCycleOffset()
 759 {
 760   // Generate randomly the strobe
 761   // phase w.r.t to the LHC clock
 762   fReadOutCycleOffset = fReadOutCycleLength*gRandom->Rndm();
 763   // fReadOutCycleOffset = 25e-9*gRandom->Rndm(); // clm: I think this way we shift too much 10-30 us! The global shift should be between the BCs?!
 764   // RS: 25 ns is too small number, the staggering will not work. Let's at the moment keep fully random shift (still, no particle from correct
 765   // collision will be lost) untill real number is specified
 766  //
 767 }
 768
 769 //______________________________________________________________________
 770 void AliITSUSimulationPix::SetResponseParam(AliITSUParamList* resp)
 771 {
 772   // attach response parameterisation data
 773   fResponseParam = resp;
 774   //
 775   int spreadID = Nint(fResponseParam->GetParameter(AliITSUSimulationPix::kChargeSpreadType));
 776   const char* hname = 0;
 777   fSpread2DHisto = 0;
 778   //
 779   switch (spreadID) {
 780     //
 781   case kSpreadFunHisto:
 782     fSpreadFun = &AliITSUSimulationPix::SpreadFrom2DHisto;
 783     hname = fResponseParam->GetParName(AliITSUSimulationPix::kChargeSpreadType);
 784     if (!(fSpread2DHisto=(TH2*)fResponseParam->GetParamObject(hname)))
 785       AliFatal(Form("Did not find 2D histo %s for charge spread parameterization",hname));
 786     break;
 787     //
 788   case kSpreadFunDoubleGauss2D:
 789     fSpreadFun = &AliITSUSimulationPix::SpreadFunDoubleGauss2D;
 790     break;
 791     //
 792   case kSpreadFunGauss2D:
 793     fSpreadFun = &AliITSUSimulationPix::SpreadFunGauss2D;
 794     break;
 795     //
 796   default: AliFatal(Form("Did not find requested spread function id=%d",spreadID));
 797   }
 798   //
 799   int readoutType = Nint(fResponseParam->GetParameter(kReadOutSchemeType));
 800   switch (readoutType) {
 801   case kReadOutStrobe:
 802     fROTimeFun = &AliITSUSimulationPix::GetReadOutCycle;
 803     break;
 804   case kReadOutRollingShutter:
 805     fROTimeFun = &AliITSUSimulationPix::GetReadOutCycleRollingShutter;
 806     break;
 807   default: AliFatal(Form("Did not find requested readout time type id=%d",readoutType));
 808   }
 809   //___ Set the Rolling Shutter read-out window
 810   fReadOutCycleLength = fResponseParam->GetParameter(kReadOutCycleLength);
 811   //___ Pixel discrimination threshold, and the S/N cut
 812   fSimuParam->SetPixThreshold(fResponseParam->GetParameter(kPixNoiseMPV) *fResponseParam->GetParameter(kPixSNDisrcCut) , fResponseParam->GetParameter(kPixSNDisrcCut),-1); //for all modules
 813   //___ Minimum number of electrons to add
 814   fSimuParam->SetPixMinElToAdd(fResponseParam->GetParameter(kPixMinElToAdd));
 815   //___ Set the Pixel Noise MPV and Sigma (the noise distribution is Landau not Gauss due to RTN)
 816   fSimuParam->SetPixNoise( fResponseParam->GetParameter(kPixNoiseMPV), fResponseParam->GetParameter(kPixNoiseSigma), -1); //for all modules
 817   //___ Pixel fake hit rate
 818   fSimuParam->SetPixFakeRate( fResponseParam->GetParameter(kPixFakeRate) );
 819   //___ To apply the noise or not
 820   if (  fResponseParam->GetParameter(kPixNoiseIsOn) > 0.01)  fSimuParam->SetPixAddNoisyFlag(kTRUE);
 821   else fSimuParam->SetPixAddNoisyFlag(kFALSE);
 822   //
 823   if(fResponseParam->GetParameter(kPixNoiseInAllMod) > 0.01 ) fSimuParam->SetPixNoiseInAllMod(kTRUE);
 824   else fSimuParam->SetPixNoiseInAllMod(kFALSE);
 825   //
 826   //  Double_t vGeVToQ = fSimuParam->GetGeVToCharge();
 827   fGlobalChargeScale = fResponseParam->GetParameter(kSpreadFunGlobalQScale);
 828
 829   AliDebug(10,Form("=============== Setting the response start ============================"));
 830   AliDebug(10,Form("=============== RO type: %d",readoutType));
 831   AliDebug(10,Form("=============== RO cycle lenght: %lf",fReadOutCycleLength));
 832   AliDebug(10,Form("=============== Noise MPV: %lf",fResponseParam->GetParameter(kPixNoiseMPV)));
 833   AliDebug(10,Form("=============== Noise Sigma: %lf",fResponseParam->GetParameter(kPixNoiseSigma)));
 834   AliDebug(10,Form("=============== Fake rate: %lf",fResponseParam->GetParameter(kPixFakeRate)));
 835   AliDebug(10,Form("=============== Noise On/Off: %d",fSimuParam->GetPixAddNoisyFlag()));
 836   AliDebug(10,Form("=============== Noise in all mod on/off: %d",fSimuParam->GetPixNoiseInAllMod()));
 837   AliDebug(10,Form("=============== Global Charge scale: %lf",fGlobalChargeScale));
 838   AliDebug(10,Form("=============== Setting the response done  ============================"));
 839
 840
 841
 842 }
 843
 844 //______________________________________________________________________
 845 Int_t AliITSUSimulationPix::GetReadOutCycleRollingShutter(Int_t row, Int_t col, Double_t hitTime)
 846 {
 847   //
 848   // Get the read-out cycle of the hit in the given column/row of the sensor.
 849   // hitTime is the time of the subhit (hit is divided to nstep charge deposit) in seconds
 850   // globalPhaseShift gives the start of the RO for the cycle in pixel wrt the LHC clock
 851   // GetRollingShutterWindow give the with of the rolling shutter read out window
 852   //
 853   double tmin = fReadOutCycleOffset + fReadOutCycleLength*(double(row)/fSeg->Npx()-1.);
 854   int cycle = Nint( (hitTime-tmin)/fReadOutCycleLength - 0.5 );
 855   AliDebug(3,Form("Rolling shutter at row%d/col%d: particle time: %e, tmin: %e : tmax %e -> cycle:%d",row,col,hitTime,tmin,
 856                   tmin+fReadOutCycleLength,cycle));
 857   return cycle;
 858   //
 859 }
 860
 861 //______________________________________________________________________
 862 Int_t AliITSUSimulationPix::GetReadOutCycle(Int_t row, Int_t col, Double_t hitTime)
 863 {
 864   //
 865   // Check whether the hit is in the read out window of the given column/row of the sensor
 866   //
 867   AliDebug(3,Form("Strobe readout: row%d/col%d: particle time: %e, tmin: %e, tmax %e",
 868                   row,col,hitTime,fReadOutCycleOffset,fReadOutCycleOffset+fReadOutCycleLength));
 869   hitTime -= fReadOutCycleOffset+0.5*fReadOutCycleLength;
 870   return (hitTime<0 || hitTime>fReadOutCycleLength) ? kMaxROCycleAccept+1 : 0;
 871   //
 872 }
 873
 874 //_______________________________________________________________________
 875 void AliITSUSimulationPix::CalcDiodeShiftInPixel(Int_t xlin, Int_t zcol, Float_t &x, Float_t &)
 876 {
 877   //
 878   // Calculates the shift of the diode wrt the geometric center of the pixel.
 879   // It is needed for staggerred pixel layout or double diode pixels with assymetric center
 880   // The shift can depend on the column or line or both...
 881   // The x and z are passed in cm
 882   //
 883
 884   TString parTitle = fResponseParam->GetTitle();
 885
 886   // M32terP31 is staggered the diode shift within pixel depends on the column
 887   if ( parTitle.Contains("M32terP31") )
 888   {
 889     if ( zcol%2 == 0 )    x += 0.30 * fSeg->Dpx(xlin);
 890     else                  x -= 0.19 * fSeg->Dpx(xlin);
 891   }
 892
 893
 894 }
 895 //_______________________________________________________________________