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