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 "AliITSUModule.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::kDetTypePix);
  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::kDetTypePix);
  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::SDigitiseModule()
 158 {
 159   //  This function begins the work of creating S-Digits.
 160
 161   AliDebug(10,Form("In event %d module %d there are %d hits", fEvent, fModule->GetIndex(),fModule->GetNHits()));
 162   if (fModule->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::FinishSDigitiseModule()
 185 {
 186    //  This function calls SDigitsToDigits which creates Digits from SDigits
 187   FrompListToDigits();
 188   ClearMap();
 189   return;
 190 }
 191
 192 //______________________________________________________________________
 193 void AliITSUSimulationPix::DigitiseModule()
 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   FinishSDigitiseModule();
 202 }
 203
 204 //______________________________________________________________________
 205 void AliITSUSimulationPix::Hits2SDigits()
 206 {
 207   // Does the charge distributions using Gaussian diffusion charge charing.
 208   Int_t nhits = fModule->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 (!fModule->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  = 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  = 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   //    AliITSUModule *mod  Pointer to this module
 283   //
 284   TObjArray *hits = fModule->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 (!fModule->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  = 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 module (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 module 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 module (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(fModule->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 = fModule->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::kDetTypePix, &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 = fModule->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   printf("added %d\n",ncand);
 626   return ncand;
 627 }
 628
 629
 630 //______________________________________________________________________
 631 void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old)
 632 {
 633   //  Take into account the coupling between adiacent pixels.
 634   //  The parameters probcol and probrow are the probability of the
 635   //  signal in one pixel shared in the two adjacent pixels along
 636   //  the column and row direction, respectively.
 637   //  Note pList is goten via GetMap() and module is not need any more.
 638   //  Otherwise it is identical to that coded by Tiziano Virgili (BSN).
 639   UInt_t col,row;
 640   Int_t iCycle;
 641   Double_t pulse1,pulse2;
 642   Double_t couplR=0.0,couplC=0.0;
 643   Double_t xr=0.;
 644   //
 645   fSensMap->GetMapIndex(old->GetUniqueID(),col,row,iCycle);
 646   int cycle = iCycle;
 647   fSimuParam->GetPixCouplingParam(couplC,couplR);
 648   if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,couplC=%e couplR=%e",
 649                                 col,row,couplC,couplR));
 650   pulse2 = pulse1 = old->GetSignal();
 651   if (pulse1<fSimuParam->GetPixMinElToAdd()) return; // too small signal
 652   for (Int_t isign=-1;isign<=1;isign+=2) {
 653     //
 654     // loop in col direction
 655     int j1 = int(col) + isign;
 656     xr = gRandom->Rndm();
 657     if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,old->GetTrack(0),old->GetHit(0),pulse1,cycle);
 658     //
 659     // loop in row direction
 660     int j2 = int(row) + isign;
 661     xr = gRandom->Rndm();
 662     if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),old->GetTrack(0),old->GetHit(0),pulse2,cycle);
 663   }
 664   //
 665 }
 666
 667 //______________________________________________________________________
 668 void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old)
 669 {
 670   //  Take into account the coupling between adiacent pixels.
 671   //  The parameters probcol and probrow are the fractions of the
 672   //  signal in one pixel shared in the two adjacent pixels along
 673   //  the column and row direction, respectively.
 674   // Inputs:
 675   // old            existing AliITSUSDigit
 676   // ntrack         track incex number
 677   // idhit          hit index number
 678   // module         module number
 679   //
 680   UInt_t col,row;
 681   Int_t cycle;
 682   Int_t modId = fModule->GetIndex();
 683   Double_t pulse1,pulse2;
 684   Double_t couplR=0.0,couplC=0.0;
 685   //
 686   fSensMap->GetMapIndex(old->GetUniqueID(),col,row,cycle);
 687   fSimuParam->GetPixCouplingParam(couplC,couplR);
 688   if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,roCycle=%d)  couplC=%e couplR=%e",col,row,cycle,couplC,couplR));
 689  //
 690  if (old->GetSignal()<fSimuParam->GetPixMinElToAdd()) return; // too small signal
 691  for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
 692    pulse2 = pulse1 = old->GetSignal();
 693    //
 694    int j1 = int(col)+isign;
 695    pulse1 *= couplC;
 696    if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(modId))) pulse1 = old->GetSignal();
 697    else UpdateMapSignal(UInt_t(j1),row,old->GetTrack(0),old->GetHit(0),pulse1,cycle);
 698
 699    // loop in row direction
 700    int j2 = int(row) + isign;
 701    pulse2 *= couplR;
 702    if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(modId))) pulse2 = old->GetSignal();
 703    else UpdateMapSignal(col,UInt_t(j2),old->GetTrack(0),old->GetHit(0),pulse2,cycle);
 704  } // for isign
 705 }
 706
 707 //______________________________________________________________________
 708 void AliITSUSimulationPix::SetResponseParam(AliITSUParamList* resp)
 709 {
 710   // attach response parameterisation data
 711   fResponseParam = resp;
 712   //
 713   int spreadID = Nint(fResponseParam->GetParameter(AliITSUSimulationPix::kChargeSpreadType));
 714   const char* hname = 0;
 715   fSpread2DHisto = 0;
 716   //
 717   switch (spreadID) {
 718     //
 719   case kSpreadFunHisto:
 720     fSpreadFun = &AliITSUSimulationPix::SpreadFrom2DHisto;
 721     hname = fResponseParam->GetParName(AliITSUSimulationPix::kChargeSpreadType);
 722     if (!(fSpread2DHisto=(TH2*)fResponseParam->GetParamObject(hname)))
 723       AliFatal(Form("Did not find 2D histo %s for charge spread parameterization",hname));
 724     break;
 725     //
 726   case kSpreadFunDoubleGauss2D:
 727     fSpreadFun = &AliITSUSimulationPix::SpreadFunDoubleGauss2D;
 728     break;
 729     //
 730   case kSpreadFunGauss2D:
 731     fSpreadFun = &AliITSUSimulationPix::SpreadFunGauss2D;
 732     break;
 733     //
 734   default: AliFatal(Form("Did not find requested spread function id=%d",spreadID));
 735   }
 736   //
 737   int readoutType = Nint(fResponseParam->GetParameter(kReadOutSchemeType));
 738   switch (readoutType) {
 739   case kReadOutStrobe:
 740     fROTimeFun = &AliITSUSimulationPix::GetReadOutCycle;
 741     break;
 742   case kReadOutRollingShutter:
 743     fROTimeFun = &AliITSUSimulationPix::GetReadOutCycleRollingShutter;
 744     break;
 745   default: AliFatal(Form("Did not find requested readout time type id=%d",readoutType));
 746   }
 747   //___ Set the Rolling Shutter read-out window
 748   fReadOutCycleLength = fResponseParam->GetParameter(kReadOutCycleLength);
 749   //___ Pixel discrimination threshold, and the S/N cut
 750   fSimuParam->SetPixThreshold(fResponseParam->GetParameter(kPixNoiseMPV) *fResponseParam->GetParameter(kPixSNDisrcCut) , fResponseParam->GetParameter(kPixSNDisrcCut),-1); //for all modules
 751   //___ Minimum number of electrons to add
 752   fSimuParam->SetPixMinElToAdd(fResponseParam->GetParameter(kPixMinElToAdd));
 753   //___ Set the Pixel Noise MPV and Sigma (the noise distribution is Landau not Gauss due to RTN)
 754   fSimuParam->SetPixNoise( fResponseParam->GetParameter(kPixNoiseMPV), fResponseParam->GetParameter(kPixNoiseSigma), -1); //for all modules
 755   //___ Pixel fake hit rate
 756   fSimuParam->SetPixFakeRate( fResponseParam->GetParameter(kPixFakeRate) );
 757   //___ To apply the noise or not
 758   if (  fResponseParam->GetParameter(kPixNoiseIsOn) > 0.01)  fSimuParam->SetPixAddNoisyFlag(kTRUE);
 759   else fSimuParam->SetPixAddNoisyFlag(kFALSE);
 760   //
 761   if(fResponseParam->GetParameter(kPixNoiseInAllMod) > 0.01 ) fSimuParam->SetPixNoiseInAllMod(kTRUE);
 762   else fSimuParam->SetPixNoiseInAllMod(kFALSE);
 763   //
 764   //  Double_t vGeVToQ = fSimuParam->GetGeVToCharge();
 765   fGlobalChargeScale = fResponseParam->GetParameter(kSpreadFunGlobalQScale);
 766
 767   AliDebug(10,Form("=============== Setting the response start ============================"));
 768   AliDebug(10,Form("=============== RO type: %d",readoutType));
 769   AliDebug(10,Form("=============== RO cycle lenght: %lf",fReadOutCycleLength));
 770   AliDebug(10,Form("=============== Noise MPV: %lf",fResponseParam->GetParameter(kPixNoiseMPV)));
 771   AliDebug(10,Form("=============== Noise Sigma: %lf",fResponseParam->GetParameter(kPixNoiseSigma)));
 772   AliDebug(10,Form("=============== Fake rate: %lf",fResponseParam->GetParameter(kPixFakeRate)));
 773   AliDebug(10,Form("=============== Noise On/Off: %d",fSimuParam->GetPixAddNoisyFlag()));
 774   AliDebug(10,Form("=============== Noise in all mod on/off: %d",fSimuParam->GetPixNoiseInAllMod()));
 775   AliDebug(10,Form("=============== Global Charge scale: %lf",fGlobalChargeScale));
 776   AliDebug(10,Form("=============== Setting the response done  ============================"));
 777
 778
 779
 780 }
 781
 782 //______________________________________________________________________
 783 Int_t AliITSUSimulationPix::GetReadOutCycleRollingShutter(Int_t row, Int_t col, Double_t hitTime)
 784 {
 785   //
 786   // Get the read-out cycle of the hit in the given column/row of the sensor.
 787   // hitTime is the time of the subhit (hit is divided to nstep charge deposit) in seconds
 788   // globalPhaseShift gives the start of the RO for the cycle in pixel wrt the LHC clock
 789   // GetRollingShutterWindow give the with of the rolling shutter read out window
 790   //
 791   double tmin = fReadOutCycleOffset + fReadOutCycleLength*(double(row)/fSeg->Npx()-1.);
 792   int cycle = Nint( (hitTime-tmin)/fReadOutCycleLength - 0.5 );
 793   AliDebug(3,Form("Rolling shutter at row%d/col%d: particle time: %e, tmin: %e : tmax %e -> cycle:%d",row,col,hitTime,tmin,
 794                   tmin+fReadOutCycleLength,cycle));
 795   return cycle;
 796   //
 797 }
 798
 799 //______________________________________________________________________
 800 Int_t AliITSUSimulationPix::GetReadOutCycle(Int_t row, Int_t col, Double_t hitTime)
 801 {
 802   //
 803   // Check whether the hit is in the read out window of the given column/row of the sensor
 804   //
 805   AliDebug(3,Form("Strobe readout: row%d/col%d: particle time: %e, tmin: %e, tmax %e",
 806                   row,col,hitTime,fReadOutCycleOffset,fReadOutCycleOffset+fReadOutCycleLength));
 807   hitTime -= fReadOutCycleOffset+0.5*fReadOutCycleLength;
 808   return (hitTime<0 || hitTime>fReadOutCycleLength) ? kMaxROCycleAccept+1 : 0;
 809   //
 810 }
 811
 812 //_______________________________________________________________________
 813 void AliITSUSimulationPix::CalcDiodeShiftInPixel(Int_t xrow, Int_t zcol, Float_t &x, Float_t &z)
 814 {
 815   //
 816   // Calculates the shift of the diode wrt the geometric center of the pixel.
 817   // It is needed for staggerred pixel layout or double diode pixels with assymetric center
 818   // The shift can depend on the column or line or both...
 819   // The x and z are passed in cm
 820   //
 821   ((AliITSUSegmentationPix*)fSeg)->GetDiodShift(xrow,zcol,x,z);
 822   //
 823 }
 824 //_______________________________________________________________________