]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/AliITSUSimulationPix.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Little coverity fix
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUSimulationPix.cxx
CommitLineData
451f5018 1
2/**************************************************************************
f6b66dab 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 **************************************************************************/
451f5018 16
451f5018 17#include <TGeoGlobalMagField.h>
18#include <TH1.h>
29ad4146 19#include <TH2.h>
451f5018 20#include <TString.h>
21#include "AliITSU.h"
22#include "AliITSUDigitPix.h"
23#include "AliITSUHit.h"
852af72e 24#include "AliITSUChip.h"
451f5018 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"
29ad4146 35#include "AliITSUParamList.h"
451f5018 36
e61afb80 37using std::cout;
38using std::endl;
a11ef2e4 39using namespace TMath;
e61afb80 40
451f5018 41ClassImp(AliITSUSimulationPix)
42////////////////////////////////////////////////////////////////////////
43// Version: 1
f6b66dab 44// Modified by D. Elia, G.E. Bruno, H. Tydesjo
451f5018 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//
29ad4146 57// 2013 Feb: Added MonoPix response and nois calculation al la MIMOSA32 (levente.molnar@cern.ch)
58//
451f5018 59////////////////////////////////////////////////////////////////////////
60
61//______________________________________________________________________
62AliITSUSimulationPix::AliITSUSimulationPix()
63: fTanLorAng(0)
f6b66dab 64,fGlobalChargeScale(1.0)
65,fSpread2DHisto(0)
66,fSpreadFun(0)
67,fROTimeFun(0)
451f5018 68{
f6b66dab 69 // Default constructor.
70 SetUniqueID(AliITSUGeomTGeo::kChipTypePix);
451f5018 71}
72
73//______________________________________________________________________
74AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap* map)
f6b66dab 75:AliITSUSimulation(sim,map)
76,fTanLorAng(0)
77,fGlobalChargeScale(1.0)
78,fSpread2DHisto(0)
79,fSpreadFun(0)
80,fROTimeFun(0)
451f5018 81{
f6b66dab 82 // standard constructor
83 SetUniqueID(AliITSUGeomTGeo::kChipTypePix);
84 Init();
451f5018 85}
86
87//______________________________________________________________________
f6b66dab 88AliITSUSimulationPix::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)
451f5018 95{
f6b66dab 96 // Copy Constructor
451f5018 97}
98
99
100//______________________________________________________________________
101AliITSUSimulationPix::~AliITSUSimulationPix()
102{
f6b66dab 103 // destructor
104 // only the sens map is owned and it is deleted by ~AliITSUSimulation
451f5018 105}
106
107//______________________________________________________________________
108AliITSUSimulationPix& AliITSUSimulationPix::operator=(const AliITSUSimulationPix &s)
109{
f6b66dab 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;
451f5018 120}
121
122//______________________________________________________________________
123void AliITSUSimulationPix::Init()
124{
f6b66dab 125 // Initilization
126 if (fSimuParam->GetPixLorentzDrift()) SetTanLorAngle(fSimuParam->GetPixLorentzHoleWeight());
451f5018 127}
128
129//______________________________________________________________________
f6b66dab 130Bool_t AliITSUSimulationPix::SetTanLorAngle(Double_t weightHole)
451f5018 131{
f6b66dab 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;
451f5018 154}
155
156//_____________________________________________________________________
852af72e 157void AliITSUSimulationPix::SDigitiseChip()
451f5018 158{
f6b66dab 159 // This function begins the work of creating S-Digits.
29ad4146 160
f6b66dab 161 AliDebug(10,Form("In event %d chip %d there are %d hits", fEvent, fChip->GetIndex(),fChip->GetNHits()));
162 if (fChip->GetNHits()) {
163 if(fResponseParam->GetParameter(kDigitalSim) == 0 ) Hits2SDigitsFast(); // analogue chip response simulation
164 else Hits2SDigitsFastDigital(); // digital chip response
165 }
166 if (!fSensMap->GetEntries()) return;
167 WriteSDigits();
168 ClearMap();
169 //
451f5018 170}
171
172//______________________________________________________________________
173void AliITSUSimulationPix::WriteSDigits()
174{
f6b66dab 175 // This function adds each S-Digit to pList
176 static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
177 int nsd = fSensMap->GetEntries();
178
179
180 for (int i=0;i<nsd;i++) {
181 AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
182 if (!sd->GetSumSignal()>0 || fSensMap->IsDisabled(sd)) continue;
183 aliITS->AddSumDigit(*sd);
184 }
185 return;
451f5018 186}
187
188//______________________________________________________________________
852af72e 189void AliITSUSimulationPix::FinishSDigitiseChip()
451f5018 190{
f6b66dab 191 // This function calls SDigitsToDigits which creates Digits from SDigits
192 FrompListToDigits();
193 ClearMap();
194 return;
451f5018 195}
196
197//______________________________________________________________________
852af72e 198void AliITSUSimulationPix::DigitiseChip()
451f5018 199{
f6b66dab 200 // This function creates Digits straight from the hits and then adds
201 // electronic noise to the digits before adding them to pList
202 // Each of the input variables is passed along to Hits2SDigits
203 //
204 if(fResponseParam->GetParameter(kDigitalSim) == 0 ) Hits2SDigitsFast(); // analogue chip response simulation
205 else Hits2SDigitsFastDigital(); // digital chip response
206 FinishSDigitiseChip();
451f5018 207}
208
209//______________________________________________________________________
c92b1537 210void AliITSUSimulationPix::Hits2SDigits()
451f5018 211{
f6b66dab 212 // Does the charge distributions using Gaussian diffusion charge charing.
213 Int_t nhits = fChip->GetNHits();
214 if (!nhits) return;
215 //
216 Int_t h,ix,iz,i;
217 Int_t idtrack;
218 Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
219 Double_t tof,x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0;
220 Double_t t,tp,st,dt=0.2,el;
221 Double_t thick = 0.5*fSeg->Dy(); // Half Thickness
451f5018 222
f6b66dab 223 //
224 for (h=0;h<nhits;h++) {
225 //
226 if (!fChip->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
227 st = Sqrt(x1*x1+y1*y1+z1*z1);
228 if (st>0.0) {
229 st = (Double_t)((Int_t)(st*1e4)); // number of microns
230 if (st<=1.0) st = 1.0;
231 dt = 1.0/st; // RS TODO: do we need 1 micron steps?
232 double dy = dt*thick;
233 y = -0.5*dy;
234 for (t=0.0;t<1.0;t+=dt) { // Integrate over t
235 tp = t+0.5*dt;
236 x = x0+x1*tp;
237 y += dy;
238 z = z0+z1*tp;
239 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
240 el = fGlobalChargeScale * dt * de / fSimuParam->GetGeVToCharge();
241 //
242 if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
243 SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
244 } // end for t
245 } else { // st == 0.0 deposit it at this point
246 x = x0;
247 y = y0 + 0.5*thick;
248 z = z0;
249 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
250 el = fGlobalChargeScale * de / fSimuParam->GetGeVToCharge();
251 if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
252 SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
253 } // end if st>0.0
254 } // Loop over all hits h
255 //
256 // Coupling
257 int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
258 AliITSUSDigit* dg = 0;
259 switch (fSimuParam->GetPixCouplingOption()) {
260 case AliITSUSimuParam::kNoCouplingPix :
261 break;
262 case AliITSUSimuParam::kNewCouplingPix :
263 for (i=nd;i--;) {
264 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
265 if (fSensMap->IsDisabled(dg)) continue;
266 SetCoupling(dg);
267 }
268 break;
269 case AliITSUSimuParam::kOldCouplingPix:
270 for (i=nd;i--;) {
271 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
272 if (fSensMap->IsDisabled(dg)) continue;
273 SetCouplingOld(dg);
274 }
275 break;
276 default:
277 break;
278
279 } // end switch
280 if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
451f5018 281}
282
283//______________________________________________________________________
c92b1537 284void AliITSUSimulationPix::Hits2SDigitsFast()
451f5018 285{
f6b66dab 286 // Does the charge distributions using Gaussian diffusion charge charing. // Inputs:
287 // AliITSUChip *mod Pointer to this chip
288 //
289 TObjArray *hits = fChip->GetHits();
290 Int_t nhits = hits->GetEntriesFast();
291 if (nhits<=0) return;
292 //
293 Int_t h,ix,iz,i;
294 Int_t idtrack;
295 Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
296 Double_t tof,x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0;
297 Double_t step,st,el,de=0.0;
298 Double_t minDim = Min(fSeg->Dpx(1),fSeg->Dpz(1)); // RStmp: smallest pitch
299 Double_t thick = fSeg->Dy();
300 //
301 for (h=0;h<nhits;h++) {
302 //
303 if (!fChip->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
304 //
305 st = Sqrt(x1*x1+y1*y1+z1*z1);
306 if (st>0.0) {
307 int np = int(1.5*st/minDim); //RStmp: inject the points in such a way that there is ~1.5 point per cell
308 np = TMath::Max(1.0*np,fResponseParam->GetParameter(kSpreadFunMinSteps));
309 AliDebug(10,Form(" Number of charge injection steps is set to %d ",np));
310 double dstep = 1./np;
311 double dy = dstep*thick;
312 y = -0.5*dy;
313 step = -0.5*dstep;
314 for (i=0;i<np;i++) { //RStmp Integrate over t
315 // for (i=0;i<kn10;i++) { // Integrate over t
316 step += dstep; // RStmp kti[i];
317 x = x0+x1*step;
318 y += dy;
319 z = z0+z1*step;
320 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
321 el = fGlobalChargeScale * dstep * de/fSimuParam->GetGeVToCharge();
322 if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
323 SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
324 } // end for i // End Integrate over t
325 }
326 else { // st == 0.0 deposit it at this point
327 x = x0;
328 y = y0+0.5*thick;
329 z = z0;
330 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
331 el = fGlobalChargeScale * de / fSimuParam->GetGeVToCharge();
332 if (fSimuParam->GetPixLorentzDrift()) x += y*fTanLorAng;
333 SpreadCharge2D(x,z,y,ix,iz,el,tof,idtrack,h);
334 } // end if st>0.0
335
336 } // Loop over all hits h
451f5018 337
f6b66dab 338 // Coupling
339 int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
340 AliITSUSDigit* dg = 0;
341 switch (fSimuParam->GetPixCouplingOption()) {
342 case AliITSUSimuParam::kNoCouplingPix :
343 break;
344 case AliITSUSimuParam::kNewCouplingPix :
345 for (i=nd;i--;) {
346 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
347 if (fSensMap->IsDisabled(dg)) continue;
348 SetCoupling(dg);
349 }
350 case AliITSUSimuParam::kOldCouplingPix:
351 for (i=nd;i--;) {
352 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
353 if (fSensMap->IsDisabled(dg)) continue;
354 SetCouplingOld(dg);
355 }
356 break;
357 default:
358 break;
359 } // end switch
360 if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
451f5018 361}
362
363//______________________________________________________________________
c92b1537 364void AliITSUSimulationPix::SpreadCharge2D(Double_t x0,Double_t z0, Double_t dy, Int_t ix0,Int_t iz0,
f6b66dab 365 Double_t el, Double_t tof, Int_t tID, Int_t hID)
451f5018 366{
f6b66dab 367 // Spreads the charge over neighboring cells. Assume charge is distributed
368 // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
369 // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
370 // Defined this way, the integral over all x and z is el.
371 // Inputs:
372 // Double_t x0 x position of point where charge is liberated (local)
373 // Double_t z0 z position of point where charge is liberated (local)
374 // Double_t dy distance from the entrance surface (diffusion sigma may depend on it)
375 // Int_t ix0 row of cell corresponding to point x0
376 // Int_t iz0 columb of cell corresponding to point z0
377 // Double_t el number of electrons liberated in this step
378 // Double_t sigx Sigma difusion along x for this step (y0 dependent)
379 // Double_t sigz Sigma difusion along z for this step (z0 dependent)
380 // Int_t tID track number
381 // Int_t hID hit "hit" index number
382 //
383 Int_t ix,iz,ixs,ixe,izs,ize;
384 Float_t x,z; // keep coordinates float (required by AliSegmentation)
385 Float_t xdioshift = 0 , zdioshift = 0 ;
386 Double_t s,dtIn[kNDtSpread]; // data transfered to spread function for integral calculation
387 //
388 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));
389 //
390 Double_t &x1 = dtIn[kCellX1];
391 Double_t &x2 = dtIn[kCellX2];
392 Double_t &z1 = dtIn[kCellZ1];
393 Double_t &z2 = dtIn[kCellZ2];
394 //
395 int nx = GetResponseParam()->GetParameter(kSpreadFunParamNXoffs);
396 int nz = GetResponseParam()->GetParameter(kSpreadFunParamNZoffs);
397 //
398 dtIn[kCellYDepth] = dy;
399 ixs = Max(-nx+ix0,0);
400 ixe = Min( nx+ix0,fSeg->Npx()-1);
401 izs = Max(-nz+iz0,0);
402 ize = Min( nz+iz0,fSeg->Npz()-1);
403 for (ix=ixs;ix<=ixe;ix++)
404 for (iz=izs;iz<=ize;iz++) {
405 //
406 // Check if the hit is inside readout window
407 int cycleRO = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof);
408 if (Abs(cycleRO)>kMaxROCycleAccept) continue;
409 //
410 fSeg->DetToLocal(ix,iz,x,z); // pixel center
411 xdioshift = zdioshift = 0;
412 double dxi = fSeg->Dpx(ix);
413 double dzi = fSeg->Dpz(iz);
414 CalcDiodeShiftInPixel(ix,iz,xdioshift,zdioshift); // Check and apply diode shift if needed
415 xdioshift *= dxi;
416 zdioshift *= dzi;
417 dxi *= 0.5;
418 dzi *= 0.5;
419 // printf("DShift: %d %d -> %.4f %.4f\n",ix,iz,xdioshift,zdioshift);
420 x1 = (x + xdioshift) - x0; // calculate distance of cell boundaries from injection center
421 z1 = (z + zdioshift) - z0;
422 x2 = x1 + dxi; // Upper
423 x1 -= dxi; // Lower
424 z2 = z1 + dzi; // Upper
425 z1 -= dzi; // Lower
426 s = el* (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fSpreadFun)(dtIn);
427 if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,tID,hID,s,cycleRO);
428 } // end for ix, iz
429 //
451f5018 430}
431
c92b1537 432//______________________________________________________________________
433Double_t AliITSUSimulationPix::SpreadFunDoubleGauss2D(const Double_t *dtIn)
434{
f6b66dab 435 // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
436 // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
437 // The spread function is assumed to be double gaussian in 2D
438 // Parameters should be: mean0,sigma0, mean1,sigma1, relative strenght of 2nd gaussian wrt 1st one
439 //
440 // 1st gaussian
441 double intg1 = GausInt2D(fResponseParam->GetParameter(kG2SigX0), // sigX
442 dtIn[kCellX1]-fResponseParam->GetParameter(kG2MeanX0), // x1-xmean
443 dtIn[kCellX2]-fResponseParam->GetParameter(kG2MeanX0), // x2-xmean
444 fResponseParam->GetParameter(kG2SigZ0), // sigZ
445 dtIn[kCellZ1]-fResponseParam->GetParameter(kG2MeanZ0), // z1-zmean
446 dtIn[kCellZ2]-fResponseParam->GetParameter(kG2MeanZ0)); // z2-zmean
447 // 2nd gaussian
448 double intg2 = GausInt2D(fResponseParam->GetParameter(kG2SigX1), // sigX
449 dtIn[kCellX1]-fResponseParam->GetParameter(kG2MeanX1), // x1-xmean
450 dtIn[kCellX2]-fResponseParam->GetParameter(kG2MeanX1), // x2-xmean
451 fResponseParam->GetParameter(kG2SigZ1), // sigZ
452 dtIn[kCellZ1]-fResponseParam->GetParameter(kG2MeanZ1), // z1-zmean
453 dtIn[kCellZ2]-fResponseParam->GetParameter(kG2MeanZ1)); // z2-zmean
454 double scl = fResponseParam->GetParameter(kG2ScaleG2);
455 return (intg1+intg2*scl)/(1+scl);
456 //
457}
c92b1537 458
29ad4146 459
460//______________________________________________________________________
461Double_t AliITSUSimulationPix::SpreadFrom2DHisto(const Double_t *dtIn)
462{
f6b66dab 463 // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
464 // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
465 // The spread function integral is taken from fSpread2DHisto extracted from the sensor response parameters
466 // list in the method SetResponseParam. The histo must return the fraction of charge integrates in the
467 // cell whose center is on the distance X=(dtIn[kCellX1]+dtIn[kCellX2])/2 and Z=(dtIn[kCellZ1]+dtIn[kCellZ2])/2
468 // from the injection point.
469 //
470 Double_t qpixfrac = 0;
471 Double_t xintp = 1e4*(dtIn[kCellX1]+dtIn[kCellX2])/2.0;
472 Double_t zintp = 1e4*(dtIn[kCellZ1]+dtIn[kCellZ2])/2.0;
473 //
474 qpixfrac = fSpread2DHisto->Interpolate(xintp,zintp); //the PSF map is given in um but the dtIn is in cm so we need to convert it
475 //
476 return qpixfrac;
29ad4146 477}
478
c92b1537 479//______________________________________________________________________
480Double_t AliITSUSimulationPix::SpreadFunGauss2D(const Double_t *dtIn)
481{
f6b66dab 482 // calculate integral of charge in the cell with boundaries at X=dtIn[kCellX1]:dtIn[kCellX2]
483 // and Z=dtIn[kCellZ1]:dtIn[kCellZ2]
484 // The spread function is assumed to be gaussian in 2D
485 // Parameters should be: mean0,sigma0
486 return GausInt2D(fResponseParam->GetParameter(kG1SigX), // sigX
487 dtIn[kCellX1]-fResponseParam->GetParameter(kG1MeanX), // x1-xmean
488 dtIn[kCellX2]-fResponseParam->GetParameter(kG1MeanX), // x2-xmean
489 fResponseParam->GetParameter(kG1SigZ), // sigZ
490 dtIn[kCellZ1]-fResponseParam->GetParameter(kG1MeanZ), // z1-zmean
491 dtIn[kCellZ2]-fResponseParam->GetParameter(kG1MeanZ)); // z2-zmean
492 //
493}
c92b1537 494
451f5018 495//______________________________________________________________________
f6b66dab 496void AliITSUSimulationPix::RemoveDeadPixels()
451f5018 497{
f6b66dab 498 // Removes dead pixels on each chip (ladder)
499 // This should be called before going from sdigits to digits (i.e. from FrompListToDigits)
500
501 AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibDead();
502 if (!calObj) return;
503 //
504 if (calObj->IsBad()) {ClearMap(); return;} // whole chip is masked
505 //
506 // prepare the list of r/o cycles seen
507 Char_t cyclesSeen[2*kMaxROCycleAccept+1];
508 int ncycles = 0;
509 for (int i=(2*kMaxROCycleAccept+1);i--;) if (fCyclesID[i]) cyclesSeen[ncycles++]=i-kMaxROCycleAccept;
510
511 // remove single bad pixels one by one
512 int nsingle = calObj->GetNrBadSingle();
513 UInt_t col,row;
514 Int_t cycle;
515 for (int i=nsingle;i--;) {
516 calObj->GetBadPixelSingle(i,row,col);
517 for (int icl=ncycles;icl--;) fSensMap->DeleteItem(col,row,cyclesSeen[icl]);
518 }
519 int nsd = fSensMap->GetEntriesUnsorted();
520 for (int isd=nsd;isd--;) {
521 AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->AtUnsorted(isd);
522 if (fSensMap->IsDisabled(sd)) continue;
523 fSensMap->GetMapIndex(sd->GetUniqueID(),col,row,cycle);
524 int chip = fSeg->GetChipFromChannel(0,col);
525 // if (calObj->IsChipMarkedBad(chip)) fSensMap->Disable(sd); // this will simple mark the hit as bad
526 if (calObj->IsChipMarkedBad(chip)) fSensMap->DeleteItem(sd); // this will suppress hit in the sorted list
527 }
528 //
451f5018 529}
530
531//______________________________________________________________________
f6b66dab 532void AliITSUSimulationPix::AddNoisyPixels()
451f5018 533{
f6b66dab 534 // Adds noisy pixels on each chip (ladder)
535 // This should be called before going from sdigits to digits (i.e. FrompListToDigits)
536 AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibNoisy();
537 if (!calObj) { AliDebug(10,Form(" No Calib Object for Noise!!! ")); return;}
538 for (Int_t i=calObj->GetNrBad(); i--;)
539 {
540 if ( fResponseParam->GetParameter(kDigitalSim) < 1.0 )
541 UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i),10*fSimuParam->GetPixThreshold(fChip->GetIndex()));
542 else
543 UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i),kNoisyPixOCDB );
544 }
545 //
451f5018 546}
547
548//______________________________________________________________________
f6b66dab 549void AliITSUSimulationPix::FrompListToDigits()
451f5018 550{
f6b66dab 551 // add noise and electronics, perform the zero suppression and add the digits to the list
552 //
553 // RS may use for ordered random sample generation dl.acm.org/ft_gateway.cfm?id=356313&type=pdf
554 //
555 int nsd = fSensMap->GetEntriesUnsorted(); // sdigits added from the signal
556 //
557 // add different kinds of noise.
558 Bool_t addNoisy = fSimuParam->GetPixAddNoisyFlag() && (nsd>0 || fSimuParam->GetPixNoiseInAllMod()); // do we generate noise?
559 if (addNoisy) {
560 AddNoisyPixels(); // constantly noisy channels
561 AddRandomNoisePixels(0.0); // random noise: at the moment generate noise only for instance 0
562 nsd = fSensMap->GetEntriesUnsorted();
451f5018 563 }
f6b66dab 564 //
565 if (nsd && fSimuParam->GetPixRemoveDeadFlag()) {
566 RemoveDeadPixels();
567 // note that here we shall use GetEntries instead of GetEntriesUnsorted since the
568 // later operates on the array where the elements are not removed by flagged
569 nsd = fSensMap->GetEntries();
451f5018 570 }
f6b66dab 571 if (!nsd) return; // nothing to digitize
572 //
573 UInt_t row,col;
574 Int_t iCycle,modId = fChip->GetIndex();
575 Double_t sig;
576 const Int_t knmaxtrk=AliITSdigit::GetNTracks();
577 static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
578 static AliITSUDigitPix dig;
579 //
580 for (int i=0;i<nsd;i++) {
581 AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
582 if (fSensMap->IsDisabled(sd)) continue;
583 //
69c15b90 584 sig=sd->GetSumSignal();
f6b66dab 585 if ( fResponseParam->GetParameter(kDigitalSim) < 1.0 &&
7413ca5b 586 sig<=fSimuParam->GetPixThreshold(modId)) continue; //Threshold only applies in analogue simulation
f6b66dab 587 //
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(),col,row,iCycle);
594 dig.SetCoord1(col);
595 dig.SetCoord2(row);
596 dig.SetROCycle(iCycle);
597 dig.SetSignal((Int_t)sig);
598 dig.SetSignalPix((Int_t)sig);
599 int ntr = sd->GetNTracks();
600 for (int j=0;j<ntr;j++) {
601 dig.SetTrack(j,sd->GetTrack(j));
602 dig.SetHit(j,sd->GetHit(j));
603 }
604 for (int j=ntr;j<knmaxtrk;j++) {
605 dig.SetTrack(j,-3);
606 dig.SetHit(j,-1);
607 }
608 aliITS->AddSimDigit(AliITSUGeomTGeo::kChipTypePix, &dig);
451f5018 609 }
f6b66dab 610 //
451f5018 611}
612
613//______________________________________________________________________
f6b66dab 614Int_t AliITSUSimulationPix::AddRandomNoisePixels(Double_t tof)
b2679935 615{
d8b69f72 616 // create random noisy sdigits above threshold
617 //
618 int modId = fChip->GetIndex();
619 int npix = fSeg->GetNPads();
620 int ncand = gRandom->Poisson( npix*fSimuParam->GetPixFakeRate() );
621 if (ncand<1) return 0;
622 //
623 double probNoisy,noiseSig,noiseMean,thresh;
624 //
625 UInt_t row,col;
626 Int_t iCycle;
627 static TArrayI ordSampleInd(100),ordSample(100); //RS!!! static is not thread-safe!!!
628 ncand = GenOrderedSample(npix,ncand,ordSample,ordSampleInd);
629 int* ordV = ordSample.GetArray();
630 int* ordI = ordSampleInd.GetArray();
631 //
632 if ( fResponseParam->GetParameter(kDigitalSim) < 1.0 ) {
633 thresh = fSimuParam->GetPixThreshold(modId);
634 fSimuParam->GetPixNoise(modId, noiseSig, noiseMean);
635 probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold
636 //
637 for (int j=0;j<ncand;j++) {
638 fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],col,row,iCycle); // create noisy digit
639 iCycle = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(row,col,tof);
640 UpdateMapNoise(col,row,AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,noiseMean,noiseSig), iCycle);
f6b66dab 641 }
d8b69f72 642 }
643 else {
644 for (int j=0;j<ncand;j++) {
645 fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],col,row,iCycle); // create noisy digit
646 iCycle = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(row,col,tof);
647 UpdateMapNoise(col,row,kNoisyPixRnd, iCycle);
f6b66dab 648 }
d8b69f72 649 }
650 return ncand;
451f5018 651}
652
1cdd0456 653
451f5018 654//______________________________________________________________________
f6b66dab 655void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old)
451f5018 656{
f6b66dab 657 // Take into account the coupling between adiacent pixels.
658 // The parameters probcol and probrow are the probability of the
659 // signal in one pixel shared in the two adjacent pixels along
660 // the column and row direction, respectively.
661 // Note pList is goten via GetMap() and chip is not need any more.
662 // Otherwise it is identical to that coded by Tiziano Virgili (BSN).
663 UInt_t col,row;
664 Int_t iCycle;
665 Double_t pulse1,pulse2;
666 Double_t couplR=0.0,couplC=0.0;
667 Double_t xr=0.;
668 //
669 fSensMap->GetMapIndex(old->GetUniqueID(),col,row,iCycle);
670 int cycle = iCycle;
671 fSimuParam->GetPixCouplingParam(couplC,couplR);
672 if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,couplC=%e couplR=%e",
673 col,row,couplC,couplR));
674 pulse2 = pulse1 = old->GetSignal();
675 if (pulse1<fSimuParam->GetPixMinElToAdd()) return; // too small signal
676 for (Int_t isign=-1;isign<=1;isign+=2) {
677 //
678 // loop in col direction
679 int j1 = int(col) + isign;
680 xr = gRandom->Rndm();
681 if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,old->GetTrack(0),old->GetHit(0),pulse1,cycle);
682 //
683 // loop in row direction
684 int j2 = int(row) + isign;
685 xr = gRandom->Rndm();
686 if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),old->GetTrack(0),old->GetHit(0),pulse2,cycle);
687 }
688 //
451f5018 689}
690
691//______________________________________________________________________
f6b66dab 692void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old)
451f5018 693{
f6b66dab 694 // Take into account the coupling between adiacent pixels.
695 // The parameters probcol and probrow are the fractions of the
696 // signal in one pixel shared in the two adjacent pixels along
697 // the column and row direction, respectively.
698 // Inputs:
699 // old existing AliITSUSDigit
700 // ntrack track incex number
701 // idhit hit index number
702 // chip chip number
703 //
704 UInt_t col,row;
705 Int_t cycle;
706 Int_t modId = fChip->GetIndex();
707 Double_t pulse1,pulse2;
708 Double_t couplR=0.0,couplC=0.0;
709 //
710 fSensMap->GetMapIndex(old->GetUniqueID(),col,row,cycle);
711 fSimuParam->GetPixCouplingParam(couplC,couplR);
712 if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,roCycle=%d) couplC=%e couplR=%e",col,row,cycle,couplC,couplR));
713 //
714 if (old->GetSignal()<fSimuParam->GetPixMinElToAdd()) return; // too small signal
715 for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
716 pulse2 = pulse1 = old->GetSignal();
717 //
718 int j1 = int(col)+isign;
719 pulse1 *= couplC;
720 if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(modId))) pulse1 = old->GetSignal();
721 else UpdateMapSignal(UInt_t(j1),row,old->GetTrack(0),old->GetHit(0),pulse1,cycle);
722
723 // loop in row direction
724 int j2 = int(row) + isign;
725 pulse2 *= couplR;
726 if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(modId))) pulse2 = old->GetSignal();
727 else UpdateMapSignal(col,UInt_t(j2),old->GetTrack(0),old->GetHit(0),pulse2,cycle);
728 } // for isign
451f5018 729}
730
c92b1537 731//______________________________________________________________________
29ad4146 732void AliITSUSimulationPix::SetResponseParam(AliITSUParamList* resp)
c92b1537 733{
f6b66dab 734 // attach response parameterisation data
735 fResponseParam = resp;
736 //
737 int spreadID = Nint(fResponseParam->GetParameter(AliITSUSimulationPix::kChargeSpreadType));
738 const char* hname = 0;
739 fSpread2DHisto = 0;
740 //
741 switch (spreadID) {
742 //
743 case kSpreadFunHisto:
744 fSpreadFun = &AliITSUSimulationPix::SpreadFrom2DHisto;
745 hname = fResponseParam->GetParName(AliITSUSimulationPix::kChargeSpreadType);
746 if (!(fSpread2DHisto=(TH2*)fResponseParam->GetParamObject(hname)))
747 AliFatal(Form("Did not find 2D histo %s for charge spread parameterization",hname));
748 break;
749 //
750 case kSpreadFunDoubleGauss2D:
751 fSpreadFun = &AliITSUSimulationPix::SpreadFunDoubleGauss2D;
752 break;
753 //
754 case kSpreadFunGauss2D:
755 fSpreadFun = &AliITSUSimulationPix::SpreadFunGauss2D;
756 break;
757 //
758 default: AliFatal(Form("Did not find requested spread function id=%d",spreadID));
759 }
760 //
761 int readoutType = Nint(fResponseParam->GetParameter(kReadOutSchemeType));
762 switch (readoutType) {
763 case kReadOutStrobe:
764 fROTimeFun = &AliITSUSimulationPix::GetReadOutCycle;
765 break;
766 case kReadOutRollingShutter:
767 fROTimeFun = &AliITSUSimulationPix::GetReadOutCycleRollingShutter;
768 break;
769 default: AliFatal(Form("Did not find requested readout time type id=%d",readoutType));
770 }
771
772 //___ Set the Rolling Shutter read-out window
773 fReadOutCycleLength = fResponseParam->GetParameter(kReadOutCycleLength);
774 //___ Pixel discrimination threshold, and the S/N cut
775 fSimuParam->SetPixThreshold(fResponseParam->GetParameter(kPixNoiseMPV) *fResponseParam->GetParameter(kPixSNDisrcCut) , fResponseParam->GetParameter(kPixSNDisrcCut),-1); //for all chips
776 //___ Minimum number of electrons to add
777 fSimuParam->SetPixMinElToAdd(fResponseParam->GetParameter(kPixMinElToAdd));
778 //___ Set the Pixel Noise MPV and Sigma (the noise distribution is Landau not Gauss due to RTN)
779 fSimuParam->SetPixNoise( fResponseParam->GetParameter(kPixNoiseMPV), fResponseParam->GetParameter(kPixNoiseSigma), -1); //for all chips
780 //___ Pixel fake hit rate
781 fSimuParam->SetPixFakeRate( fResponseParam->GetParameter(kPixFakeRate) );
782 //___ To apply the noise or not
783 if ( fResponseParam->GetParameter(kPixNoiseIsOn) > 0.01) fSimuParam->SetPixAddNoisyFlag(kTRUE);
784 else fSimuParam->SetPixAddNoisyFlag(kFALSE);
785 //
786 if(fResponseParam->GetParameter(kPixNoiseInAllMod) > 0.01 ) fSimuParam->SetPixNoiseInAllMod(kTRUE);
787 else fSimuParam->SetPixNoiseInAllMod(kFALSE);
788 //
789 // Double_t vGeVToQ = fSimuParam->GetGeVToCharge();
790 fGlobalChargeScale = fResponseParam->GetParameter(kSpreadFunGlobalQScale);
791
792 AliDebug(10,Form("=============== Setting the response start ============================"));
793 AliDebug(10,Form("=============== Digital (1) / Analogue (0) simu: %f",fResponseParam->GetParameter(kDigitalSim)));
794 AliDebug(10,Form("=============== RO type: %d",readoutType));
795 AliDebug(10,Form("=============== RO cycle lenght: %lf",fReadOutCycleLength));
796 AliDebug(10,Form("=============== Noise MPV: %lf",fResponseParam->GetParameter(kPixNoiseMPV)));
797 AliDebug(10,Form("=============== Noise Sigma: %lf",fResponseParam->GetParameter(kPixNoiseSigma)));
798 AliDebug(10,Form("=============== Fake rate: %lf",fResponseParam->GetParameter(kPixFakeRate)));
799 AliDebug(10,Form("=============== Noise On/Off: %d",fSimuParam->GetPixAddNoisyFlag()));
800 AliDebug(10,Form("=============== Noise in all mod on/off: %d",fSimuParam->GetPixNoiseInAllMod()));
801 AliDebug(10,Form("=============== Global Charge scale: %lf",fGlobalChargeScale));
802 AliDebug(10,Form("=============== Setting the response done ============================"));
29ad4146 803
c92b1537 804}
0ebc85cf 805
806//______________________________________________________________________
b2679935 807Int_t AliITSUSimulationPix::GetReadOutCycleRollingShutter(Int_t row, Int_t col, Double_t hitTime)
0ebc85cf 808{
f6b66dab 809 //
810 // Get the read-out cycle of the hit in the given column/row of the sensor.
811 // hitTime is the time of the subhit (hit is divided to nstep charge deposit) in seconds
812 // globalPhaseShift gives the start of the RO for the cycle in pixel wrt the LHC clock
813 // GetRollingShutterWindow give the with of the rolling shutter read out window
814 //
815 double tmin = fReadOutCycleOffset + fReadOutCycleLength*(double(row)/fSeg->Npx()-1.);
816 int cycle = Nint( (hitTime-tmin)/fReadOutCycleLength - 0.5 );
817 AliDebug(3,Form("Rolling shutter at row%d/col%d: particle time: %e, tmin: %e : tmax %e -> cycle:%d",row,col,hitTime,tmin,
818 tmin+fReadOutCycleLength,cycle));
819 return cycle;
820 //
0ebc85cf 821}
822
823//______________________________________________________________________
b2679935 824Int_t AliITSUSimulationPix::GetReadOutCycle(Int_t row, Int_t col, Double_t hitTime)
0ebc85cf 825{
f6b66dab 826 //
827 // Check whether the hit is in the read out window of the given column/row of the sensor
828 //
829 AliDebug(3,Form("Strobe readout: row%d/col%d: particle time: %e, tmin: %e, tmax %e",
830 row,col,hitTime,fReadOutCycleOffset,fReadOutCycleOffset+fReadOutCycleLength));
831 hitTime -= fReadOutCycleOffset+0.5*fReadOutCycleLength;
832 return (hitTime<0 || hitTime>fReadOutCycleLength) ? kMaxROCycleAccept+1 : 0;
833 //
0ebc85cf 834}
835
29ad4146 836//_______________________________________________________________________
ee58ce21 837void AliITSUSimulationPix::CalcDiodeShiftInPixel(Int_t xrow, Int_t zcol, Float_t &x, Float_t &z)
29ad4146 838{
f6b66dab 839 //
840 // Calculates the shift of the diode wrt the geometric center of the pixel.
841 // It is needed for staggerred pixel layout or double diode pixels with assymetric center
842 // The shift can depend on the column or line or both...
843 // The x and z are passed in cm
844 //
845 ((AliITSUSegmentationPix*)fSeg)->GetDiodShift(xrow,zcol,x,z);
846 //
29ad4146 847}
f6b66dab 848
849//______________________________________________________________________
850void AliITSUSimulationPix::Hits2SDigitsFastDigital()
851{
852 // Does the digital chip response simulation
853 // AliITSUChip *mod Pointer to this chip
854 //
855
856
857 TObjArray *hits = fChip->GetHits();
858 Int_t nhits = hits->GetEntriesFast();
859 if (nhits<=0) return;
860 //
861 Int_t h,ix,iz;
862 Int_t idtrack;
f6b66dab 863 Double_t tof,x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0;
69c15b90 864 Double_t el,de=0.0;
f6b66dab 865
866 //
867 for (h=0;h<nhits;h++) {
868 //
869 if (!fChip->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,tof,idtrack)) continue;
870 //
69c15b90 871 //double st = Sqrt(x1*x1+y1*y1+z1*z1);
f6b66dab 872
873 //___ place hit to the middle of the path segment - CHANGE LATER !
69c15b90 874 // keep coordinates float (required by AliSegmentation)
875 float x = (x0+x1)/2.0;
876 //float y = (y0+y1)/2.0;
877 float z = (z0+z1)/2.0;
f6b66dab 878 //
879 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
880 //
881 // Check if the hit is inside readout window
882 int cycleRO = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof);
883 if (Abs(cycleRO)>kMaxROCycleAccept) continue;
884 //
885 //
886 el = de / fSimuParam->GetGeVToCharge();
887 //
888 PlaceDigitalPixels(x,z,el,tof,idtrack,h);
889 //
890 } // Loop over all hits h
891
892}
893//______________________________________________________________________
894void AliITSUSimulationPix::PlaceDigitalPixels(Double_t x0hit,Double_t z0hit, Double_t el, Double_t tof, Int_t tID, Int_t hID)
895{
896 // Place the digital pixel positions on the sensor
897 // Inputs:
898 // Double_t x0hit x position of point where charge is liberated (local) - hit
899 // Double_t z0hit z position of point where charge is liberated (local) - hit
900 // Double_t el number of electrons liberated in this step
901 // Double_t sigx Sigma difusion along x for this step (y0 dependent)
902 // Double_t sigz Sigma difusion along z for this step (z0 dependent)
903 // Int_t tID track number
904 // Int_t hID hit "hit" index number
905 //
906
907
908 Int_t ix,iz,nx,nz;
909 Float_t x,z; // keep coordinates float (required by AliSegmentation)
910 Float_t distX = 0, distZ = 0;
911
912 //___ TEMPORARY - place a fixed pattern cluster COG to a distance d(x,z) away from hit - TEMPORARY
913 // Pattern used (not realistic ebye but averaging over events yes):
914 // -+-
915 // +++
916 // -+-
917 //
918 //___ This list should come from a look up table based on CluTypeID as well as COG coord
919
920 //
921 TRandom3 rnd;
922 distX = rnd.Uniform(-5.0*1e-4, 5.0*1e-4); //in um
923 distZ = rnd.Uniform(-5.0*1e-4, 5.0*1e-4); //in um
924 //
925 x = x0hit + distX;
926 z = z0hit + distZ;
927 //
928 if(!fSeg->LocalToDet(x,z,ix,iz)) return; // if clu CoG is outside of the chip skipp the cluster -> refine later
929 //
930 const Int_t nCluPixels = 5;
931 Int_t aPixListX[nCluPixels] = { 0, -1, 0, 1, 0};
932 Int_t aPixListZ[nCluPixels] = { 1, 0, 0, 0, -1};
933 //
934 Double_t s = el / 1.0 / nCluPixels;
935 //
936 int cycleRO;
937 //
938 for (Int_t ipix = 0 ; ipix < nCluPixels; ipix++)
939 {
940 nx = ix + aPixListX[ipix];
941 nz = iz + aPixListZ[ipix];
942 cycleRO = (((AliITSUSimulationPix*)this)->*AliITSUSimulationPix::fROTimeFun)(ix,iz,tof);
943 if ( nx >= 0 && nx <= fSeg -> Npx() && nz >= 0 && nz <= fSeg -> Npz() ) UpdateMapSignal(nz,nx,tID,hID,s,cycleRO); //if the pixel is in the detector
944 }
945
946
947
948}
949
950
951