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 returning.", 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         // Check if the hit is inside readout window
 259         if ( !(((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof) ) continue;
 260         SpreadCharge2D(x,z,y,ix,iz,el,idtrack,h);
 261       } // end for t
 262     } else { // st == 0.0 deposit it at this point
 263       x   = x0;
 264       y   = y0 + 0.5*thick;
 265       z   = z0;
 266       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 267       el  = de / fSimuParam->GetGeVToCharge();
 268       if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
 269       // Check if the hit is inside readout window
 270       if ( !(((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof) ) continue;
 271       SpreadCharge2D(x,z,y,ix,iz,el,idtrack,h);
 272     } // end if st>0.0
 273   } // Loop over all hits h
 274   //
 275   // Coupling
 276   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
 277   AliITSUSDigit* dg = 0;
 278   switch (fSimuParam->GetPixCouplingOption()) {
 279   case AliITSUSimuParam::kNewCouplingPix :
 280     for (i=nd;i--;) {
 281       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 282       if (fSensMap->IsDisabled(dg)) continue;
 283       SetCoupling(dg,idtrack,h);
 284     }
 285     break;
 286   case AliITSUSimuParam::kOldCouplingPix:
 287     for (i=nd;i--;) {
 288       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 289       if (fSensMap->IsDisabled(dg)) continue;
 290       SetCouplingOld(dg,idtrack,h);
 291     }
 292     break;
 293   default:
 294     break;
 295
 296   } // end switch
 297   if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
 298 }
 299
 300 //______________________________________________________________________
 301 void AliITSUSimulationPix::Hits2SDigitsFast()
 302 {
 303   // Does the charge distributions using Gaussian diffusion charge charing.    // Inputs:
 304   //    AliITSUModule *mod  Pointer to this module
 305   //
 306   TObjArray *hits = fModule->GetHits();
 307   Int_t nhits = hits->GetEntriesFast();
 308   if (nhits<=0) return;
 309   //
 310   Int_t h,ix,iz,i;
 311   Int_t idtrack;
 312   Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
 313   Double_t tof,x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0;
 314   Double_t step,st,el,de=0.0;
 315   Double_t minDim = Min(fSeg->Dpx(1),fSeg->Dpz(1)); // RStmp: smallest pitch
 316   Double_t thick = fSeg->Dy();
 317   //
 318   for (h=0;h<nhits;h++) {
 319     //
 320     if (!fModule->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
 321     //
 322     st = Sqrt(x1*x1+y1*y1+z1*z1);
 323     if (st>0.0) {
 324       int np = int(1.5*st/minDim);  //RStmp: inject the points in such a way that there is ~1.5 point per cell
 325       np = TMath::Max(1.0*np,fResponseParam->GetParameter(kSpreadFunMinSteps));
 326       AliDebug(10,Form(" Number of charge injection steps is set to %d ",np));
 327       double dstep = 1./np;
 328       double dy = dstep*thick;
 329       y = -0.5*dy;
 330       step = -0.5*dstep;
 331       for (i=0;i<np;i++) {          //RStmp Integrate over t
 332         //      for (i=0;i<kn10;i++) { // Integrate over t
 333         step  += dstep;  // RStmp kti[i];
 334         x   = x0+x1*step;
 335         y  += dy;
 336         z   = z0+z1*step;
 337         if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 338         el  = fGlobalChargeScale * dstep * de/fSimuParam->GetGeVToCharge();
 339         if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
 340         // Check if the hit is inside readout window
 341         if ( !(((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof) ) continue;
 342         SpreadCharge2D(x,z,y,ix,iz,el,idtrack,h);
 343       } // end for i // End Integrate over t
 344     }
 345     else { // st == 0.0 deposit it at this point
 346       x   = x0;
 347       y   = y0+0.5*thick;
 348       z   = z0;
 349       if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
 350       el  = de / fSimuParam->GetGeVToCharge();
 351       if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
 352       // Check if the hit is inside readout window
 353       if ( !(((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof) ) continue;
 354       SpreadCharge2D(x,z,y,ix,iz,el,idtrack,h);
 355     } // end if st>0.0
 356
 357   } // Loop over all hits h
 358
 359   // Coupling
 360   int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
 361   AliITSUSDigit* dg = 0;
 362   switch (fSimuParam->GetPixCouplingOption()) {
 363   case AliITSUSimuParam::kNewCouplingPix :
 364     for (i=nd;i--;) {
 365       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 366       if (fSensMap->IsDisabled(dg)) continue;
 367       SetCoupling(dg,idtrack,h);
 368     }
 369   case AliITSUSimuParam::kOldCouplingPix:
 370     for (i=nd;i--;) {
 371       dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
 372       if (fSensMap->IsDisabled(dg)) continue;
 373       SetCouplingOld(dg,idtrack,h);
 374     }
 375     break;
 376   default:
 377     break;
 378   } // end switch
 379   if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
 380 }
 381
 382 //______________________________________________________________________
 383 void AliITSUSimulationPix::SpreadCharge2D(Double_t x0,Double_t z0, Double_t dy, Int_t ix0,Int_t iz0,
 384                                           Double_t el, Int_t tID, Int_t hID)
 385 {
 386   // Spreads the charge over neighboring cells. Assume charge is distributed
 387   // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
 388   // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
 389   // Defined this way, the integral over all x and z is el.
 390   // Inputs:
 391   //    Double_t x0   x position of point where charge is liberated (local)
 392   //    Double_t z0   z position of point where charge is liberated (local)
 393   //    Double_t dy   distance from the entrance surface (diffusion sigma may depend on it)
 394   //    Int_t    ix0  row of cell corresponding to point x0
 395   //    Int_t    iz0  columb of cell corresponding to point z0
 396   //    Double_t el   number of electrons liberated in this step
 397   //    Double_t sigx Sigma difusion along x for this step (y0 dependent)
 398   //    Double_t sigz Sigma difusion along z for this step (z0 dependent)
 399   //    Int_t    tID  track number
 400   //    Int_t    hID  hit "hit" index number
 401   //
 402   Int_t ix,iz,ixs,ixe,izs,ize;
 403   Float_t x,z;   // keep coordinates float (required by AliSegmentation)
 404   Float_t xdioshift = 0 , zdioshift = 0 ;
 405   Double_t s,dtIn[kNDtSpread]; // data transfered to spread function for integral calculation
 406   //
 407   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));
 408   //
 409   Double_t &x1 = dtIn[kCellX1];
 410   Double_t &x2 = dtIn[kCellX2];
 411   Double_t &z1 = dtIn[kCellZ1];
 412   Double_t &z2 = dtIn[kCellZ2];
 413   //
 414   int nx = GetResponseParam()->GetParameter(kSpreadFunParamNXoffs);
 415   int nz = GetResponseParam()->GetParameter(kSpreadFunParamNZoffs);
 416   //
 417   dtIn[kCellYDepth]  = dy;
 418   ixs = Max(-nx+ix0,0);
 419   ixe = Min( nx+ix0,fSeg->Npx()-1);
 420   izs = Max(-nz+iz0,0);
 421   ize = Min( nz+iz0,fSeg->Npz()-1);
 422   for (ix=ixs;ix<=ixe;ix++)
 423     for (iz=izs;iz<=ize;iz++) {
 424       fSeg->DetToLocal(ix,iz,x,z); // pixel center
 425       xdioshift = zdioshift = 0;
 426       CalcDiodeShiftInPixel(ix,iz,xdioshift,zdioshift);    // Check and apply diode shift if needed
 427       double dxi = 0.5*fSeg->Dpx(ix);
 428       double dzi = 0.5*fSeg->Dpz(iz);
 429       x1  = (x + xdioshift) - x0;   // calculate distance of cell boundaries from injection center
 430       z1  = (z + zdioshift) - z0;
 431       x2  = x1 + dxi; // Upper
 432       x1 -= dxi;      // Lower
 433       z2  = z1 + dzi; // Upper
 434       z1 -= dzi;      // Lower
 435       s   = el* (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fSpreadFun)(dtIn);
 436       if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,tID,hID,s);
 437     } // end for ix, iz
 438   //
 439 }
 440
 441 //______________________________________________________________________
 442 Double_t AliITSUSimulationPix::SpreadFunDoubleGauss2D(const Double_t *dtIn)
 443 {
 444   // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
 445   // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
 446   // The spread function is assumed to be double gaussian in 2D
 447   // Parameters should be: mean0,sigma0, mean1,sigma1, relative strenght of 2nd gaussian wrt 1st one
 448   //
 449   // 1st gaussian
 450   double intg1 = GausInt2D(fResponseParam->GetParameter(kG2SigX0),  // sigX
 451                            dtIn[kCellX1]-fResponseParam->GetParameter(kG2MeanX0),      // x1-xmean
 452                            dtIn[kCellX2]-fResponseParam->GetParameter(kG2MeanX0),      // x2-xmean
 453                            fResponseParam->GetParameter(kG2SigZ0),  // sigZ
 454                            dtIn[kCellZ1]-fResponseParam->GetParameter(kG2MeanZ0),    // z1-zmean
 455                            dtIn[kCellZ2]-fResponseParam->GetParameter(kG2MeanZ0));   // z2-zmean
 456   // 2nd gaussian
 457   double intg2 = GausInt2D(fResponseParam->GetParameter(kG2SigX1),  // sigX
 458                            dtIn[kCellX1]-fResponseParam->GetParameter(kG2MeanX1),      // x1-xmean
 459                            dtIn[kCellX2]-fResponseParam->GetParameter(kG2MeanX1),      // x2-xmean
 460                            fResponseParam->GetParameter(kG2SigZ1),  // sigZ
 461                            dtIn[kCellZ1]-fResponseParam->GetParameter(kG2MeanZ1),    // z1-zmean
 462                            dtIn[kCellZ2]-fResponseParam->GetParameter(kG2MeanZ1));   // z2-zmean
 463   double scl = fResponseParam->GetParameter(kG2ScaleG2);
 464   return (intg1+intg2*scl)/(1+scl);
 465   //
 466 }
 467
 468
 469 //______________________________________________________________________
 470 Double_t AliITSUSimulationPix::SpreadFrom2DHisto(const Double_t *dtIn)
 471 {
 472   // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
 473   // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
 474   // The spread function integral is taken from fSpread2DHisto extracted from the sensor response parameters
 475   // list in the method SetResponseParam. The histo must return the fraction of charge integrates in the
 476   // cell whose center is on the distance X=(dtIn[kCellX1]+dtIn[kCellX2])/2 and Z=(dtIn[kCellZ1]+dtIn[kCellZ2])/2
 477   // from the injection point.
 478   //
 479   Double_t qpixfrac = 0;
 480   Double_t xintp = 1e4*(dtIn[kCellX1]+dtIn[kCellX2])/2.0;
 481   Double_t zintp = 1e4*(dtIn[kCellZ1]+dtIn[kCellZ2])/2.0;
 482   //
 483   qpixfrac =  fSpread2DHisto->Interpolate(xintp,zintp); //the PSF map is given in um but the dtIn is in cm so we need to convert it
 484   //
 485   return qpixfrac;
 486 }
 487
 488 //______________________________________________________________________
 489 Double_t AliITSUSimulationPix::SpreadFunGauss2D(const Double_t *dtIn)
 490 {
 491   // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
 492   // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
 493   // The spread function is assumed to be gaussian in 2D
 494   // Parameters should be: mean0,sigma0
 495   return GausInt2D(fResponseParam->GetParameter(kG1SigX),  // sigX
 496                    fResponseParam->GetParameter(kG1SigZ),  // sigZ
 497                    dtIn[kCellX1]-fResponseParam->GetParameter(kG1MeanX),    // x1-xmean
 498                    dtIn[kCellX2]-fResponseParam->GetParameter(kG1MeanX),    // x2-xmean
 499                    dtIn[kCellZ1]-fResponseParam->GetParameter(kG1MeanZ),    // z1-zmean
 500                    dtIn[kCellZ2]-fResponseParam->GetParameter(kG1MeanZ));   // z2-zmean
 501   //
 502 }
 503
 504 //______________________________________________________________________
 505 void AliITSUSimulationPix::RemoveDeadPixels()
 506 {
 507   // Removes dead pixels on each module (ladder)
 508   // This should be called before going from sdigits to digits (FrompListToDigits)
 509
 510   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibDead();
 511   if (!calObj) return;
 512   //
 513   if (calObj->IsBad()) {ClearMap(); return;} // whole module is masked
 514   //
 515   // remove single bad pixels one by one
 516   int nsingle = calObj->GetNrBadSingle();
 517   UInt_t col,row;
 518   for (int i=nsingle;i--;) {
 519     calObj->GetBadPixelSingle(i,row,col);
 520     fSensMap->DeleteItem(col,row);
 521   }
 522   int nsd = fSensMap->GetEntriesUnsorted();
 523   for (int isd=nsd;isd--;) {
 524     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->AtUnsorted(isd);
 525     if (fSensMap->IsDisabled(sd)) continue;
 526     fSensMap->GetMapIndex(sd->GetUniqueID(),col,row);
 527     int chip = fSeg->GetChipFromChannel(0,col);
 528     //    if (calObj->IsChipMarkedBad(chip)) fSensMap->Disable(sd); // this will simple mark the hit as bad
 529     if (calObj->IsChipMarkedBad(chip)) fSensMap->DeleteItem(sd); // this will suppress hit in the sorted list
 530   }
 531   //
 532 }
 533
 534 //______________________________________________________________________
 535 void AliITSUSimulationPix::AddNoisyPixels()
 536 {
 537   // Adds noisy pixels on each module (ladder)
 538   // This should be called before going from sdigits to digits (FrompListToDigits)
 539   AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibNoisy();
 540   if (!calObj) { AliDebug(10,Form("  No Calib Object for Noise!!! ")); return;}
 541   for (Int_t i=calObj->GetNrBad(); i--;) UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i),
 542                                                         10*fSimuParam->GetPixThreshold(fModule->GetIndex()));
 543   //
 544 }
 545
 546 //______________________________________________________________________
 547 void AliITSUSimulationPix::FrompListToDigits()
 548 {
 549   // add noise and electronics, perform the zero suppression and add the
 550   // digit to the list
 551   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
 552   UInt_t ix,iz;
 553   Double_t sig;
 554   const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
 555   static AliITSUDigitPix dig;
 556   // RS: in principle:
 557   // 1) for every pixel w/o hit we have to generate a noise and activate the pixel if the noise exceeds the threshold.
 558   // 2) for every pixel with hit we should add random noise and check if the total signal exceeds the threshold
 559   // With many channels this will be too time consuming, hence I do the following
 560   // 1) Precalculate the probability that the nois alone will exceed the threshold (probnoisy).
 561   // 2) Chose randomly empty pixels according to fakes rate
 562   // 3) For pixels having a hits apply the usual noise and compare with threshold
 563   //
 564   // RS may use for ordered random sample generation dl.acm.org/ft_gateway.cfm?id=356313&type=pdf
 565   //
 566   int maxInd = fSensMap->GetMaxIndex();
 567   int modId = fModule->GetIndex();
 568   //
 569   int nsd = fSensMap->GetEntries();
 570   Int_t prevID=0,curID=0;
 571   TArrayI ordSampleInd(100),ordSample(100);
 572   //
 573   double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(modId);
 574   fSimuParam->GetPixNoise(modId, noiseSig, noiseMean);
 575   probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold
 576   //
 577   for (int i=0;i<nsd;i++) {
 578     AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
 579     if (fSensMap->IsDisabled(sd)) continue;
 580     curID = (int)sd->GetUniqueID();
 581     //
 582     if ( fSimuParam->GetPixAddNoisyFlag() ) {
 583       CreateNoisyDigits(prevID,curID,probNoisy, noiseSig, noiseMean);
 584       prevID = curID+1;
 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);
 594     dig.SetCoord1(iz);
 595     dig.SetCoord2(ix);
 596     dig.SetSignal((Int_t)sig);
 597     dig.SetSignalPix((Int_t)sig);
 598     int ntr = sd->GetNTracks();
 599     for (int j=0;j<ntr;j++) {
 600       dig.SetTrack(j,sd->GetTrack(j));
 601       dig.SetHit(j,sd->GetHit(j));
 602     }
 603     for (int j=ntr;j<knmaxtrk;j++) {
 604       dig.SetTrack(j,-3);
 605       dig.SetHit(j,-1);
 606     }
 607     aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix, &dig);
 608   }
 609   // if needed, add noisy pixels with id from last real hit to maxID
 610   if (fSimuParam->GetPixAddNoisyFlag() &&
 611       (nsd>0  || ( nsd==0 && fSimuParam->GetPixNoiseInAllMod()) )) {
 612     AliDebug(10,Form("CreateNoisyDigits2( %d ,%d) module %d",prevID,maxInd,modId));
 613     CreateNoisyDigits(prevID,maxInd,probNoisy, noiseSig, noiseMean);
 614   }
 615   //
 616 }
 617
 618 //______________________________________________________________________
 619 Int_t AliITSUSimulationPix::CreateNoisyDigits(Int_t minID,Int_t maxID,double probNoisy, double noise, double base)
 620 {//RS PAYATT2
 621   // create random noisy digits above threshold within id range [minID,maxID[
 622   // see FrompListToDigits for details
 623   //
 624   static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
 625   UInt_t ix,iz;
 626   static AliITSUDigitPix dig;
 627   static TArrayI ordSampleInd(100),ordSample(100); //RS!!! static is not thread-safe!!!
 628   const Int_t    knmaxtrk=AliITSdigit::GetNTracks();
 629   //
 630   Int_t ncand = 0;
 631   int npix = maxID-minID;
 632   if (npix<1 || (ncand=gRandom->Poisson(npix * fSimuParam->GetPixFakeRate() ))<1) return ncand; // decide how many noisy pixels will be added
 633   ncand = GenOrderedSample(npix,ncand,ordSample,ordSampleInd);
 634   int* ordV = ordSample.GetArray();
 635   int* ordI = ordSampleInd.GetArray();
 636   for (int j=0;j<ncand;j++) {
 637     fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],iz,ix);   // create noisy digit
 638     dig.SetCoord1(iz);
 639     dig.SetCoord2(ix);
 640     dig.SetSignal(1);
 641     dig.SetSignalPix((Int_t)AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,base,noise));
 642     for (int k=knmaxtrk;k--;) {
 643       dig.SetTrack(k,-3);
 644       dig.SetHit(k,-1);
 645     }
 646     aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix,&dig);
 647     AliDebug(10,Form("Add noisy pixel %d(%d/%d) Noise=%d",ordV[ordI[j]],iz,ix,dig.GetSignalPix()));
 648   }
 649   return ncand;
 650 }
 651
 652 //______________________________________________________________________
 653 void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old, Int_t ntrack, Int_t idhit)
 654 {
 655   //  Take into account the coupling between adiacent pixels.
 656   //  The parameters probcol and probrow are the probability of the
 657   //  signal in one pixel shared in the two adjacent pixels along
 658   //  the column and row direction, respectively.
 659   //  Note pList is goten via GetMap() and module is not need any more.
 660   //  Otherwise it is identical to that coded by Tiziano Virgili (BSN).
 661   //Begin_Html
 662   /*
 663     <img src="picts/ITS/barimodel_3.gif">
 664      </pre>
 665      <br clear=left>
 666      <font size=+2 color=red>
 667      <a href="mailto:tiziano.virgili@cern.ch"></a>.
 668      </font>
 669      <pre>
 670    */
 671    //End_Html
 672    // Inputs:
 673   // old                  existing AliITSUSDigit
 674   // Int_t ntrack         track incex number
 675   // Int_t idhit          hit index number
 676   UInt_t col,row;
 677   Double_t pulse1,pulse2;
 678   Double_t couplR=0.0,couplC=0.0;
 679   Double_t xr=0.;
 680   //
 681   fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
 682   fSimuParam->GetPixCouplingParam(couplC,couplR);
 683   if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
 684                                 col,row,ntrack,idhit,couplC,couplR));
 685   pulse2 = pulse1 = old->GetSignal();
 686   if (pulse1<fSimuParam->GetPixMinElToAdd()) return; // too small signal
 687   for (Int_t isign=-1;isign<=1;isign+=2) {
 688     //
 689     // loop in col direction
 690     int j1 = int(col) + isign;
 691     xr = gRandom->Rndm();
 692     if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
 693     //
 694     // loop in row direction
 695     int j2 = int(row) + isign;
 696     xr = gRandom->Rndm();
 697     if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
 698   }
 699   //
 700 }
 701
 702 //______________________________________________________________________
 703 void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, Int_t ntrack,Int_t idhit)
 704 {
 705   //  Take into account the coupling between adiacent pixels.
 706   //  The parameters probcol and probrow are the fractions of the
 707   //  signal in one pixel shared in the two adjacent pixels along
 708   //  the column and row direction, respectively.
 709   //Begin_Html
 710   /*
 711     <img src="picts/ITS/barimodel_3.gif">
 712     </pre>
 713     <br clear=left>
 714     <font size=+2 color=red>
 715     <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
 716     </font>
 717     <pre>
 718   */
 719   //End_Html
 720   // Inputs:
 721   // old            existing AliITSUSDigit
 722   // ntrack         track incex number
 723   // idhit          hit index number
 724   // module         module number
 725   //
 726   UInt_t col,row;
 727   Int_t modId = fModule->GetIndex();
 728   Double_t pulse1,pulse2;
 729   Double_t couplR=0.0,couplC=0.0;
 730   //
 731   fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
 732   fSimuParam->GetPixCouplingParam(couplC,couplR);
 733   if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d)  couplC=%e couplR=%e",
 734                                 col,row,ntrack,idhit,couplC,couplR));
 735  //
 736  if (old->GetSignal()<fSimuParam->GetPixMinElToAdd()) return; // too small signal
 737  for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
 738    pulse2 = pulse1 = old->GetSignal();
 739    //
 740    int j1 = int(col)+isign;
 741    pulse1 *= couplC;
 742    if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(modId))) pulse1 = old->GetSignal();
 743    else UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
 744
 745    // loop in row direction
 746    int j2 = int(row) + isign;
 747    pulse2 *= couplR;
 748    if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(modId))) pulse2 = old->GetSignal();
 749    else UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
 750  } // for isign
 751 }
 752
 753 //______________________________________________________________________
 754 void AliITSUSimulationPix::GenerateReadOutCycleOffset()
 755 {
 756   // Generate randomly the strobe
 757   // phase w.r.t to the LHC clock
 758   fReadOutCycleOffset = fReadOutCycleLength*gRandom->Rndm();
 759   // 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?!
 760   // 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
 761   // collision will be lost) untill real number is specified
 762  //
 763 }
 764
 765 //______________________________________________________________________
 766 void AliITSUSimulationPix::SetResponseParam(AliITSUParamList* resp)
 767 {
 768   // attach response parameterisation data
 769   fResponseParam = resp;
 770   //
 771   int spreadID = Nint(fResponseParam->GetParameter(AliITSUSimulationPix::kChargeSpreadType));
 772   const char* hname = 0;
 773   fSpread2DHisto = 0;
 774   //
 775   switch (spreadID) {
 776     //
 777   case kSpreadFunHisto:
 778     fSpreadFun = &AliITSUSimulationPix::SpreadFrom2DHisto;
 779     hname = fResponseParam->GetParName(AliITSUSimulationPix::kChargeSpreadType);
 780     if (!(fSpread2DHisto=(TH2*)fResponseParam->GetParamObject(hname)))
 781       AliFatal(Form("Did not find 2D histo %s for charge spread parameterization",hname));
 782     break;
 783     //
 784   case kSpreadFunDoubleGauss2D:
 785     fSpreadFun = &AliITSUSimulationPix::SpreadFunDoubleGauss2D;
 786     break;
 787     //
 788   case kSpreadFunGauss2D:
 789     fSpreadFun = &AliITSUSimulationPix::SpreadFunGauss2D;
 790     break;
 791     //
 792   default: AliFatal(Form("Did not find requested spread function id=%d",spreadID));
 793   }
 794   //
 795   int readoutType = Nint(fResponseParam->GetParameter(kReadOutSchemeType));
 796   switch (readoutType) {
 797   case kReadOutStrobe:
 798     fROTimeFun = &AliITSUSimulationPix::IsHitInReadOutWindow;
 799     break;
 800   case kReadOutRollingShutter:
 801     fROTimeFun = &AliITSUSimulationPix::IsHitInReadOutWindowRollingShutter;
 802     break;
 803   default: AliFatal(Form("Did not find requested readout time type id=%d",readoutType));
 804   }
 805   //___ Set the Rolling Shutter read-out window
 806   fReadOutCycleLength = fResponseParam->GetParameter(kReadOutCycleLength);
 807   //___ Pixel discrimination threshold, and the S/N cut
 808   fSimuParam->SetPixThreshold(fResponseParam->GetParameter(kPixNoiseMPV) *fResponseParam->GetParameter(kPixSNDisrcCut) , fResponseParam->GetParameter(kPixSNDisrcCut),-1); //for all modules
 809   //___ Minimum number of electrons to add
 810   fSimuParam->SetPixMinElToAdd(fResponseParam->GetParameter(kPixMinElToAdd));
 811   //___ Set the Pixel Noise MPV and Sigma (the noise distribution is Landau not Gauss due to RTN)
 812   fSimuParam->SetPixNoise( fResponseParam->GetParameter(kPixNoiseMPV), fResponseParam->GetParameter(kPixNoiseSigma), -1); //for all modules
 813   //___ Pixel fake hit rate
 814   fSimuParam->SetPixFakeRate( fResponseParam->GetParameter(kPixFakeRate) );
 815   //___ To apply the noise or not
 816   if (  fResponseParam->GetParameter(kPixNoiseIsOn) > 0.01)  fSimuParam->SetPixAddNoisyFlag(kTRUE);
 817   else fSimuParam->SetPixAddNoisyFlag(kFALSE);
 818   //
 819   if(fResponseParam->GetParameter(kPixNoiseInAllMod) > 0.01 ) fSimuParam->SetPixNoiseInAllMod(kTRUE);
 820   else fSimuParam->SetPixNoiseInAllMod(kFALSE);
 821   //
 822   //  Double_t vGeVToQ = fSimuParam->GetGeVToCharge();
 823   fGlobalChargeScale = fResponseParam->GetParameter(kSpreadFunGlobalQScale);
 824
 825   AliDebug(10,Form("=============== Setting the response start ============================"));
 826   AliDebug(10,Form("=============== RO type: %d",readoutType));
 827   AliDebug(10,Form("=============== RO cycle lenght: %lf",fReadOutCycleLength));
 828   AliDebug(10,Form("=============== Noise MPV: %lf",fResponseParam->GetParameter(kPixNoiseMPV)));
 829   AliDebug(10,Form("=============== Noise Sigma: %lf",fResponseParam->GetParameter(kPixNoiseSigma)));
 830   AliDebug(10,Form("=============== Fake rate: %lf",fResponseParam->GetParameter(kPixFakeRate)));
 831   AliDebug(10,Form("=============== Noise On/Off: %d",fSimuParam->GetPixAddNoisyFlag()));
 832   AliDebug(10,Form("=============== Noise in all mod on/off: %d",fSimuParam->GetPixNoiseInAllMod()));
 833   AliDebug(10,Form("=============== Global Charge scale: %lf",fGlobalChargeScale));
 834   AliDebug(10,Form("=============== Setting the response done  ============================"));
 835
 836
 837
 838 }
 839
 840 //______________________________________________________________________
 841 Bool_t AliITSUSimulationPix::IsHitInReadOutWindowRollingShutter(Int_t row, Int_t col, Double_t hitTime)
 842 {
 843   //
 844   // Check whether the hit is in the read out window of the given column/row of the sensor
 845   // hitTime is the time of the subhit (hit is divided to nstep charge deposit) in seconds
 846   // globalPhaseShift gives the start of the RO for the cycle in pixel wrt the LHC clock
 847   // GetRollingShutterWindow give the with of the rolling shutter read out window
 848   //
 849   double timePerRow = fReadOutCycleLength / fSeg->Npx();
 850   double tmax = fReadOutCycleOffset + timePerRow*(row+1);
 851   double tmin = tmax - fReadOutCycleLength;
 852   AliDebug(3,Form("Rolling shutter at row%d/col%d: particle time: %e, tmin: %e : %e",row,col,hitTime,tmin,tmax));
 853   return (hitTime<tmin || hitTime>tmax) ? kFALSE : kTRUE;
 854   //
 855 }
 856
 857 //______________________________________________________________________
 858 Bool_t AliITSUSimulationPix::IsHitInReadOutWindow(Int_t row, Int_t col, Double_t hitTime)
 859 {
 860   //
 861   // Check whether the hit is in the read out window of the given column/row of the sensor
 862   //
 863   AliDebug(3,Form("Strobe readout: row%d/col%d: particle time: %e, tmin: %e, tmax %e",
 864                   row,col,hitTime,fReadOutCycleOffset,fReadOutCycleOffset+fReadOutCycleLength));
 865   hitTime -= fReadOutCycleOffset+0.5*fReadOutCycleLength;
 866   return (hitTime<0 || hitTime>fReadOutCycleLength) ? kFALSE : kTRUE;
 867   //
 868 }
 869
 870 //_______________________________________________________________________
 871 void AliITSUSimulationPix::CalcDiodeShiftInPixel(Int_t xlin, Int_t zcol, Float_t &x, Float_t &)
 872 {
 873   //
 874   // Calculates the shift of the diode wrt the geometric center of the pixel.
 875   // It is needed for staggerred pixel layout or double diode pixels with assymetric center
 876   // The shift can depend on the column or line or both...
 877   // The x and z are passed in cm
 878   //
 879
 880   TString parTitle = fResponseParam->GetTitle();
 881
 882   // M32terP31 is staggered the diode shift within pixel depends on the column
 883   if ( parTitle.Contains("M32terP31") )
 884   {
 885     if ( zcol%2 == 0 )    x += 0.30 * fSeg->Dpx(xlin);
 886     else                  x -= 0.19 * fSeg->Dpx(xlin);
 887   }
 888
 889
 890 }
 891 //_______________________________________________________________________