ITS/UPGRADE/AliITSUSimulationPix.cxx

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