git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
suppressed std:... implicit calls
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUSimulationPix.cxx
CommitLineData
451f5018 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
451f5018 22#include <TGeoGlobalMagField.h>
23#include <TH1.h>
24#include <TString.h>
25#include "AliITSU.h"
26#include "AliITSUDigitPix.h"
27#include "AliITSUHit.h"
28#include "AliITSUModule.h"
29#include "AliITSUSensMap.h"
30#include "AliITSUCalibrationPix.h"
31#include "AliITSUSegmentationPix.h"
32#include "AliITSUSimulationPix.h"
33#include "AliLog.h"
34#include "AliRun.h"
35#include "AliMagF.h"
36#include "AliMathBase.h"
37#include "AliITSUSimuParam.h"
38#include "AliITSUSDigit.h"
39
40ClassImp(AliITSUSimulationPix)
41////////////////////////////////////////////////////////////////////////
42// Version: 1
43// Modified by D. Elia, G.E. Bruno, H. Tydesjo
44// Fast diffusion code by Bjorn S. Nilsen
45// March-April 2006
46// October 2007: GetCalibrationObjects() removed
47//
48// Version: 0
49// Written by Boris Batyunya
50// December 20 1999
51//
52// Adapted for pixels of ITS upgrade July 2012, ruben.shahoyan@cern.ch
53//
54// AliITSUSimulationPix is to do the simulation of pixels
55//
56////////////////////////////////////////////////////////////////////////
57
58//______________________________________________________________________
59AliITSUSimulationPix::AliITSUSimulationPix()
60: fTanLorAng(0)
61 ,fStrobe(kTRUE)
62 ,fStrobeLenght(4)
63 ,fStrobePhase(-12.5e-9)
64{
65 // Default constructor.
66}
67
68//______________________________________________________________________
69AliITSUSimulationPix::AliITSUSimulationPix(AliITSUSimuParam* sim,AliITSUSensMap* map)
70 :AliITSUSimulation(sim,map)
71 ,fTanLorAng(0)
72 ,fStrobe(kTRUE)
73 ,fStrobeLenght(4)
74 ,fStrobePhase(-12.5e-9)
75{
76 // standard constructor
77 Init();
78}
79
80//______________________________________________________________________
81AliITSUSimulationPix::AliITSUSimulationPix(const AliITSUSimulationPix &s)
82 :AliITSUSimulation(s)
83 ,fTanLorAng(s.fTanLorAng)
84 ,fStrobe(s.fStrobe)
85 ,fStrobeLenght(s.fStrobeLenght)
86 ,fStrobePhase(s.fStrobePhase)
87{
88 // Copy Constructor
89}
90
91
92//______________________________________________________________________
93AliITSUSimulationPix::~AliITSUSimulationPix()
94{
95 // destructor
96 // only the sens map is owned and it is deleted by ~AliITSUSimulation
97}
98
99//______________________________________________________________________
100AliITSUSimulationPix& AliITSUSimulationPix::operator=(const AliITSUSimulationPix &s)
101{
102 // Assignment operator
103 if (&s == this) return *this;
104 AliITSUSimulation::operator=(s);
105 fStrobe = s.fStrobe;
106 fStrobeLenght = s.fStrobeLenght;
107 fStrobePhase = s.fStrobePhase;
108 return *this;
109}
110
111//______________________________________________________________________
112void AliITSUSimulationPix::Init()
113{
114 // Initilization
4fa9d550 115 if (fSimuParam->GetPixLorentzDrift()) SetTanLorAngle(fSimuParam->GetPixLorentzHoleWeight());
451f5018 116}
117
118//______________________________________________________________________
119Bool_t AliITSUSimulationPix::SetTanLorAngle(Double_t weightHole)
120{
121 // This function set the Tangent of the Lorentz angle.
122 // A weighted average is used for electrons and holes
123 // Input: Double_t weightHole: wheight for hole: it should be in the range [0,1]
124 // output: Bool_t : kTRUE in case of success
125 //
126 if (weightHole<0) {
127 weightHole=0.;
128 AliWarning("You have asked for negative Hole weight");
129 AliWarning("I'm going to use only electrons");
130 }
131 if (weightHole>1) {
132 weightHole=1.;
133 AliWarning("You have asked for weight > 1");
134 AliWarning("I'm going to use only holes");
135 }
136 Double_t weightEle=1.-weightHole;
137 AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
138 if (!fld) AliFatal("The field is not initialized");
139 Double_t bz = fld->SolenoidField();
140 fTanLorAng = TMath::Tan(weightHole*fSimuParam->LorentzAngleHole(bz) +
141 weightEle*fSimuParam->LorentzAngleElectron(bz));
142 return kTRUE;
143}
144
145//_____________________________________________________________________
146void AliITSUSimulationPix::SDigitiseModule(AliITSUModule *mod, Int_t /*mask*/,Int_t event, AliITSsegmentation* seg)
147{
148 // This function begins the work of creating S-Digits.
149 if (!(mod->GetNHits())) {
150 AliDebug(1,Form("In event %d module %d there are %d hits returning.",
151 event, mod->GetIndex(),mod->GetNHits()));
152 return;
153 }
154 //
155 if (fStrobe) if (event != fEvent) GenerateStrobePhase();
156 InitSimulationModule(mod->GetIndex(), event, seg );
157 // Hits2SDigits(mod);
158 Hits2SDigitsFast(mod);
4fa9d550 159 if (fSimuParam->GetPixAddNoisyFlag()) AddNoisyPixels();
160 if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();
451f5018 161 WriteSDigits();
162 ClearMap();
163}
164
165//______________________________________________________________________
166void AliITSUSimulationPix::WriteSDigits()
167{
168 // This function adds each S-Digit to pList
169 static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
170 int nsd = fSensMap->GetEntries();
171 for (int i=0;i<nsd;i++) {
172 AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
173 if (!sd->GetSumSignal()>0 || fSensMap->IsDisabled(sd)) continue;
174 aliITS->AddSumDigit(*sd);
175 }
176 return;
177}
178
179//______________________________________________________________________
180void AliITSUSimulationPix::FinishSDigitiseModule()
181{
182 // This function calls SDigitsToDigits which creates Digits from SDigits
183 FrompListToDigits();
184 ClearMap();
185 return;
186}
187
188//______________________________________________________________________
189void AliITSUSimulationPix::DigitiseModule(AliITSUModule *mod,Int_t /*mask*/, Int_t event, AliITSsegmentation* seg)
190{
191 // This function creates Digits straight from the hits and then adds
192 // electronic noise to the digits before adding them to pList
193 // Each of the input variables is passed along to Hits2SDigits
194 //
195 if (fStrobe) if (event != fEvent) GenerateStrobePhase();
196 InitSimulationModule( mod->GetIndex(), event, seg );
197 // Hits2SDigits(mod);
198 Hits2SDigitsFast(mod);
199 //
4fa9d550 200 if (fSimuParam->GetPixAddNoisyFlag()) AddNoisyPixels();
201 if (fSimuParam->GetPixRemoveDeadFlag()) RemoveDeadPixels();
451f5018 202 FrompListToDigits();
203 ClearMap();
204}
205
206//______________________________________________________________________
207void AliITSUSimulationPix::Hits2SDigits(AliITSUModule *mod)
208{
209 // Does the charge distributions using Gaussian diffusion charge charing.
210 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
211 const Double_t kcmtomic = 1.e4;
212 const Double_t kBunchLenght = 25e-9; // LHC clock
213 Int_t nhits = mod->GetNHits();
214 if (!nhits) return;
215 //
216 Int_t h,ix,iz,i;
217 Int_t idtrack;
4fa9d550 218 Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
451f5018 219 Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0,ld=0.0;
4fa9d550 220 Double_t t,tp,st,dt=0.2,el,sig,sigx,sigz,fda;
451f5018 221 Double_t thick = 0.5*kmictocm*fSeg->Dy(); // Half Thickness
4fa9d550 222 fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
451f5018 223 //
224 for (h=0;h<nhits;h++) {
225 if (fStrobe &&
226 ((mod->GetHit(h)->GetTOF()<fStrobePhase) ||
227 (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
228 ) continue;
229 //
230 if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
231 st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
232 if (st>0.0) {
233 st = (Double_t)((Int_t)(st*kcmtomic)); // number of microns
234 if (st<=1.0) st = 1.0;
235 dt = 1.0/st;
236 for (t=0.0;t<1.0;t+=dt) { // Integrate over t
237 tp = t+0.5*dt;
238 x = x0+x1*tp;
239 y = y0+y1*tp;
240 z = z0+z1*tp;
241 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
242 el = dt * de / fSimuParam->GetGeVToCharge();
243 //
451f5018 244 sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
245 sigx=sig;
246 sigz=sig*fda;
4fa9d550 247 if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
451f5018 248 SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
249 } // end for t
250 } else { // st == 0.0 deposit it at this point
251 x = x0;
252 y = y0;
253 z = z0;
254 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
255 el = de / fSimuParam->GetGeVToCharge();
256 sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
257 sigx = sig;
258 sigz = sig*fda;
4fa9d550 259 if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
451f5018 260 SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
261 } // end if st>0.0
262 } // Loop over all hits h
263 //
264 // Coupling
265 int nd = fSensMap->GetEntriesUnsorted(); // use unsorted access when possible, since it is faster
266 AliITSUSDigit* dg = 0;
4fa9d550 267 switch (fSimuParam->GetPixCouplingOption()) {
268 case AliITSUSimuParam::kNewCouplingPix :
451f5018 269 for (i=nd;i--;) {
270 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
271 if (fSensMap->IsDisabled(dg)) continue;
272 SetCoupling(dg,idtrack,h);
273 }
274 break;
4fa9d550 275 case AliITSUSimuParam::kOldCouplingPix:
451f5018 276 for (i=nd;i--;) {
277 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
278 if (fSensMap->IsDisabled(dg)) continue;
279 SetCouplingOld(dg,idtrack,h);
280 }
281 break;
282 default:
283 break;
284
285 } // end switch
4fa9d550 286 if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
451f5018 287}
288
289//______________________________________________________________________
290void AliITSUSimulationPix::Hits2SDigitsFast(AliITSUModule *mod)
291{
292 // Does the charge distributions using Gaussian diffusion charge charing. // Inputs:
293 // AliITSUModule *mod Pointer to this module
294 // Output:
295 // none.
296 // Return:
297 // none.
298 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
299 const Int_t kn10=10;
300 const Double_t kti[kn10]={7.443716945e-3,2.166976971e-1,3.397047841e-1,
301 4.325316833e-1,4.869532643e-1,5.130467358e-1,
302 5.674683167e-1,6.602952159e-1,7.833023029e-1,
303 9.255628306e-1};
304 const Double_t kwi[kn10]={1.477621124e-1,1.346333597e-1,1.095431813e-1,
305 7.472567455e-2,3.333567215e-2,3.333567215e-2,
306 7.472567455e-2,1.095431813e-1,1.346333597e-1,
307 1.477621124e-1};
308 const Double_t kBunchLenght = 25e-9; // LHC clock
309 TObjArray *hits = mod->GetHits();
310 Int_t nhits = hits->GetEntriesFast();
311 if (nhits<=0) return;
312 //
313 Int_t h,ix,iz,i;
314 Int_t idtrack;
4fa9d550 315 Float_t x,y,z; // keep coordinates float (required by AliSegmentation)
316 Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0;
317 Double_t t,st,el,sig,sigx,sigz,fda,de=0.0,ld=0.0;
451f5018 318 Double_t thick = 0.5*kmictocm*fSeg->Dy(); // Half thickness
4fa9d550 319 fSimuParam->GetPixSigmaDiffusionAsymmetry(fda);
451f5018 320 //
321 for (h=0;h<nhits;h++) {
322 // Check if the hit is inside readout window
323 if (fStrobe &&
324 ((mod->GetHit(h)->GetTOF()<fStrobePhase) ||
325 (mod->GetHit(h)->GetTOF()>(fStrobePhase+fStrobeLenght*kBunchLenght)))
326 ) continue;
327 //
328 if (!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
329 st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
330 if (st>0.0)
331 for (i=0;i<kn10;i++) { // Integrate over t
332 t = kti[i];
333 x = x0+x1*t;
334 y = y0+y1*t;
335 z = z0+z1*t;
336 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
337 el = kwi[i]*de/fSimuParam->GetGeVToCharge();
338 sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
339 sigx=sig;
340 sigz=sig*fda;
4fa9d550 341 if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
451f5018 342 SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,idtrack,h);
451f5018 343 } // end for i // End Integrate over t
344 else { // st == 0.0 deposit it at this point
345 x = x0;
346 y = y0;
347 z = z0;
348 if (!(fSeg->LocalToDet(x,z,ix,iz))) continue; // outside
349 el = de / fSimuParam->GetGeVToCharge();
350 sig = fSimuParam->SigmaDiffusion1D(TMath::Abs(thick + y));
351 sigx=sig;
352 sigz=sig*fda;
4fa9d550 353 if (fSimuParam->GetPixLorentzDrift()) ld=(y+thick)*fTanLorAng;
451f5018 354 SpreadChargeAsym(x,z,ix,iz,el,sigx,sigz,ld,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;
4fa9d550 362 switch (fSimuParam->GetPixCouplingOption()) {
363 case AliITSUSimuParam::kNewCouplingPix :
451f5018 364 for (i=nd;i--;) {
365 dg = (AliITSUSDigit*)fSensMap->AtUnsorted(i);
366 if (fSensMap->IsDisabled(dg)) continue;
367 SetCoupling(dg,idtrack,h);
368 }
4fa9d550 369 case AliITSUSimuParam::kOldCouplingPix:
451f5018 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
4fa9d550 379 if (GetDebug(2)) AliInfo(Form("Finished fCoupling=%d",fSimuParam->GetPixCouplingOption()));
451f5018 380}
381
382//______________________________________________________________________
383void AliITSUSimulationPix::SpreadCharge(Double_t x0,Double_t z0,
384 Int_t ix0,Int_t iz0,
385 Double_t el,Double_t sig,Double_t ld,
386 Int_t t,Int_t hi)
387{
388 // Spreads the charge over neighboring cells. Assume charge is distributed
389 // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg)
390 // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig)
4fa9d550 391 // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
451f5018 392 // Defined this way, the integral over all x and z is el.
393 // Inputs:
394 // Double_t x0 x position of point where charge is liberated
395 // Double_t z0 z position of point where charge is liberated
396 // Int_t ix0 row of cell corresponding to point x0
397 // Int_t iz0 columb of cell corresponding to point z0
398 // Double_t el number of electrons liberated in this step
399 // Double_t sig Sigma difusion for this step (y0 dependent)
400 // Double_t ld lorentz drift in x for this step (y0 dependent)
401 // Int_t t track number
402 // Int_t ti hit track index number
403 // Int_t hi hit "hit" index number
404 // Outputs:
405 // none.
406 // Return:
407 // none.
408 const Int_t knx = 3,knz = 2;
409 const Double_t kRoot2 = 1.414213562; // Sqrt(2).
410 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
411 Int_t ix,iz,ixs,ixe,izs,ize;
4fa9d550 412 Float_t x,z; // keep coordinates float (required by AliSegmentation)
413 Double_t s,sp,x1,x2,z1,z2;
451f5018 414 //
4fa9d550 415 if (GetDebug(2)) AliInfo(Form("(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,sig=%e,t=%d,i=%d)",x0,z0,ix0,iz0,el,sig,t,hi));
451f5018 416 if (sig<=0.0) { // if sig<=0 No diffusion to simulate.
417 UpdateMapSignal(iz0,ix0,t,hi,el);
451f5018 418 return;
419 } // end if
420 sp = 1.0/(sig*kRoot2);
451f5018 421 ixs = TMath::Max(-knx+ix0,0);
422 ixe = TMath::Min(knx+ix0,fSeg->Npx()-1);
423 izs = TMath::Max(-knz+iz0,0);
424 ize = TMath::Min(knz+iz0,fSeg->Npz()-1);
425 for (ix=ixs;ix<=ixe;ix++)
426 for (iz=izs;iz<=ize;iz++) {
427 fSeg->DetToLocal(ix,iz,x,z); // pixel center
4fa9d550 428 double dxi = 0.5*kmictocm*fSeg->Dpx(ix);
429 double dzi = 0.5*kmictocm*fSeg->Dpz(iz);
451f5018 430 x1 = x;
431 z1 = z;
4fa9d550 432 x2 = x1 + dxi; // Upper
433 x1 -= dxi; // Lower
434 z2 = z1 + dzi; // Upper
435 z1 -= dzi; // Lower
451f5018 436 x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
437 x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
438 z1 -= z0; // Distance from where track traveled
439 z2 -= z0; // Distance from where track traveled
4fa9d550 440 s = el*0.25; // Correction based on definision of Erfc
451f5018 441 s *= AliMathBase::ErfcFast(sp*x1) - AliMathBase::ErfcFast(sp*x2);
451f5018 442 s *= AliMathBase::ErfcFast(sp*z1) - AliMathBase::ErfcFast(sp*z2);
4fa9d550 443 if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
451f5018 444 } // end for ix, iz
445 //
446}
447
448//______________________________________________________________________
449void AliITSUSimulationPix::SpreadChargeAsym(Double_t x0,Double_t z0,
4fa9d550 450 Int_t ix0,Int_t iz0,
451 Double_t el,Double_t sigx,Double_t sigz,
452 Double_t ld,Int_t t,Int_t hi)
451f5018 453{
454 // Spreads the charge over neighboring cells. Assume charge is distributed
455 // as charge(x,z) = (el/2*pi*sigx*sigz)*exp(-arg)
456 // arg=((x-x0)*(x-x0)/2*sigx*sigx)+((z-z0*z-z0)/2*sigz*sigz)
4fa9d550 457 // if fSimuParam->GetPixLorentzDrift()=kTRUE, then x0=x0+ld (Lorentz drift taken into account)
451f5018 458 // Defined this way, the integral over all x and z is el.
459 // Inputs:
460 // Double_t x0 x position of point where charge is liberated
461 // Double_t z0 z position of point where charge is liberated
462 // Int_t ix0 row of cell corresponding to point x0
463 // Int_t iz0 columb of cell corresponding to point z0
464 // Double_t el number of electrons liberated in this step
465 // Double_t sigx Sigma difusion along x for this step (y0 dependent)
4fa9d550 466 // Double_t sigz Sigma difusion along z for this step (z0 dependent)
451f5018 467 // Double_t ld lorentz drift in x for this stip (y0 dependent)
468 // Int_t t track number
469 // Int_t ti hit track index number
470 // Int_t hi hit "hit" index number
471 // Outputs:
472 // none.
473 // Return:
474 // none.
4fa9d550 475 const Int_t knx = 3,knz = 3; // RS: TO TUNE
451f5018 476 const Double_t kRoot2 = 1.414213562; // Sqrt(2).
477 const Double_t kmictocm = 1.0e-4; // convert microns to cm.
478 Int_t ix,iz,ixs,ixe,izs,ize;
4fa9d550 479 Float_t x,z; // keep coordinates float (required by AliSegmentation)
480 Double_t s,spx,spz,x1,x2,z1,z2;
451f5018 481 //
cf457606 482 if (GetDebug(2)) AliInfo(Form("(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,sigx=%e, sigz=%e, t=%d,i=%d,ld=%e)",
483 x0,z0,ix0,iz0,el,sigx,sigz,t,hi,ld));
451f5018 484 if (sigx<=0.0 || sigz<=0.0) { // if sig<=0 No diffusion to simulate.
485 UpdateMapSignal(iz0,ix0,t,hi,el);
451f5018 486 return;
487 } // end if
488 spx = 1.0/(sigx*kRoot2);
489 spz = 1.0/(sigz*kRoot2);
451f5018 490 ixs = TMath::Max(-knx+ix0,0);
491 ixe = TMath::Min(knx+ix0,fSeg->Npx()-1);
492 izs = TMath::Max(-knz+iz0,0);
493 ize = TMath::Min(knz+iz0,fSeg->Npz()-1);
494 for (ix=ixs;ix<=ixe;ix++)
495 for (iz=izs;iz<=ize;iz++) {
496 fSeg->DetToLocal(ix,iz,x,z); // pixel center
4fa9d550 497 double dxi = 0.5*kmictocm*fSeg->Dpx(ix);
498 double dzi = 0.5*kmictocm*fSeg->Dpz(iz);
451f5018 499 x1 = x;
500 z1 = z;
4fa9d550 501 x2 = x1 + dxi; // Upper
502 x1 -= dxi; // Lower
503 z2 = z1 + dzi; // Upper
504 z1 -= dzi; // Lower
451f5018 505 x1 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
506 x2 -= x0+ld; // Distance from where track traveled (taking into account the Lorentz drift)
507 z1 -= z0; // Distance from where track traveled
508 z2 -= z0; // Distance from where track traveled
4fa9d550 509 s = el*0.25; // Correction based on definision of Erfc
451f5018 510 s *= AliMathBase::ErfcFast(spx*x1) - AliMathBase::ErfcFast(spx*x2);
451f5018 511 s *= AliMathBase::ErfcFast(spz*z1) - AliMathBase::ErfcFast(spz*z2);
4fa9d550 512 if (s>fSimuParam->GetPixMinElToAdd()) UpdateMapSignal(iz,ix,t,hi,s);
451f5018 513 } // end for ix, iz
514 //
515}
516
517//______________________________________________________________________
518void AliITSUSimulationPix::RemoveDeadPixels()
519{
520 // Removes dead pixels on each module (ladder)
521 // This should be called before going from sdigits to digits (FrompListToDigits)
522
523 AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibDead();
524 if (!calObj) return;
525 //
526 if (calObj->IsBad()) {ClearMap(); return;} // whole module is masked
527 //
528 // remove single bad pixels one by one
529 int nsingle = calObj->GetNrBadSingle();
530 UInt_t col,row;
531 for (int i=nsingle;i--;) {
532 calObj->GetBadPixelSingle(i,row,col);
533 fSensMap->DeleteItem(col,row);
534 }
535 int nsd = fSensMap->GetEntriesUnsorted();
536 for (int isd=nsd;isd--;) {
537 AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->AtUnsorted(isd);
538 if (fSensMap->IsDisabled(sd)) continue;
539 fSensMap->GetMapIndex(sd->GetUniqueID(),col,row);
540 int chip = fSeg->GetChipFromChannel(0,col);
541 // if (calObj->IsChipMarkedBad(chip)) fSensMap->Disable(sd); // this will simple mark the hit as bad
542 if (calObj->IsChipMarkedBad(chip)) fSensMap->DeleteItem(sd); // this will suppress hit in the sorted list
543 }
544 //
545}
546
547//______________________________________________________________________
548void AliITSUSimulationPix::AddNoisyPixels()
549{
550 // Adds noisy pixels on each module (ladder)
551 // This should be called before going from sdigits to digits (FrompListToDigits)
552 AliITSUCalibrationPix* calObj = (AliITSUCalibrationPix*) GetCalibNoisy();
553 if (!calObj) return;
554 for (Int_t i=calObj->GetNrBad(); i--;) UpdateMapNoise(calObj->GetBadColAt(i), calObj->GetBadRowAt(i),
4fa9d550 555 10*fSimuParam->GetPixThreshold(fModule));
451f5018 556 //
557}
558
559//______________________________________________________________________
560void AliITSUSimulationPix::FrompListToDigits()
561{
562 // add noise and electronics, perform the zero suppression and add the
563 // digit to the list
564 static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
565 UInt_t ix,iz;
566 Double_t sig;
567 const Int_t knmaxtrk=AliITSdigit::GetNTracks();
568 static AliITSUDigitPix dig;
569 // RS: in principle:
570 // 1) for every pixel w/o hit we have to generate a noise and activate the pixel if the noise exceeds the threshold.
571 // 2) for every pixel with hit we should add random noise and check if the total signal exceeds the threshold
572 // With many channels this will be too time consuming, hence I do the following
573 // 1) Precalculate the probability that the nois alone will exceed the threshold.
574 // 2) Chose randomly empty pixels according to this probability and apply the noise above threshold.
575 // 3) For pixels having a hits apply the usual noise and compare with threshold
576 //
577 // RS may use for ordered random sample generation dl.acm.org/ft_gateway.cfm?id=356313&type=pdf
578 //
579 int maxInd = fSensMap->GetMaxIndex();
580 double minProb = 0.1/maxInd;
581 //
582 int nsd = fSensMap->GetEntries();
583 Int_t prevID=0,curID=0;
584 TArrayI ordSampleInd(100),ordSample(100);
585 //
4fa9d550 586 double probNoisy,noiseSig,noiseMean,thresh = fSimuParam->GetPixThreshold(fModule);
587 fSimuParam->GetPixNoise(fModule, noiseSig, noiseMean);
451f5018 588 probNoisy = AliITSUSimuParam::CalcProbNoiseOverThreshold(noiseMean,noiseSig,thresh); // prob. to have noise above threshold
589 //
590 for (int i=0;i<nsd;i++) {
591 AliITSUSDigit* sd = (AliITSUSDigit*)fSensMap->At(i); // ordered in index
592 if (fSensMap->IsDisabled(sd)) continue;
593 curID = (int)sd->GetUniqueID();
594 //
595 if (probNoisy>minProb) { // generate randomly noisy pixels above the threshold, with ID's between previous hit and current
596 CreateNoisyDigits(prevID,curID,probNoisy, noiseSig, noiseMean);
597 prevID = curID+1;
598 }
599 //
4fa9d550 600 if ((sig=sd->GetSumSignal())<=fSimuParam->GetPixThreshold(fModule)) continue;
451f5018 601 if (TMath::Abs(sig)>2147483647.0) { //RS?
602 //PH 2147483647 is the max. integer
603 //PH This apparently is a problem which needs investigation
604 AliWarning(Form("Too big or too small signal value %f",sig));
605 }
606 fSensMap->GetMapIndex(sd->GetUniqueID(),iz,ix);
607 dig.SetCoord1(iz);
608 dig.SetCoord2(ix);
609 dig.SetSignal(1);
610 dig.SetSignalPix((Int_t)sig);
611 int ntr = sd->GetNTracks();
612 for (int j=0;j<ntr;j++) {
613 dig.SetTrack(j,sd->GetTrack(j));
614 dig.SetHit(j,sd->GetHit(j));
615 }
616 for (int j=ntr;j<knmaxtrk;j++) {
617 dig.SetTrack(j,-3);
618 dig.SetHit(j,-1);
619 }
02d6eccc 620 aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix, &dig);
451f5018 621 }
622 // if needed, add noisy pixels with id from last real hit to maxID
623 if (probNoisy>minProb) CreateNoisyDigits(prevID,maxInd,probNoisy, noiseSig, noiseMean);
624 //
625}
626
627//______________________________________________________________________
628Int_t AliITSUSimulationPix::CreateNoisyDigits(Int_t minID,Int_t maxID,double probNoisy, double noise, double base)
629{
630 // create random noisy digits above threshold within id range [minID,maxID[
631 // see FrompListToDigits for details
632 //
633 static AliITSU *aliITS = (AliITSU*)gAlice->GetModule("ITS");
634 UInt_t ix,iz;
635 static AliITSUDigitPix dig;
636 static TArrayI ordSampleInd(100),ordSample(100); //RS!!! static is not thread-safe!!!
637 const Int_t knmaxtrk=AliITSdigit::GetNTracks();
638 //
639 Int_t ncand = 0;
640 int npix = maxID-minID;
641 if (npix<1 || (ncand=gRandom->Poisson(npix*probNoisy))<1) return ncand; // decide how many noisy pixels will be added
642 ncand = GenOrderedSample(npix,ncand,ordSample,ordSampleInd);
643 int* ordV = ordSample.GetArray();
644 int* ordI = ordSampleInd.GetArray();
645 for (int j=0;j<ncand;j++) {
646 fSensMap->GetMapIndex((UInt_t)ordV[ordI[j]],iz,ix); // create noisy digit
647 dig.SetCoord1(iz);
648 dig.SetCoord2(ix);
649 dig.SetSignal(1);
650 dig.SetSignalPix((Int_t)AliITSUSimuParam::GenerateNoiseQFunction(probNoisy,base,noise));
651 for (int k=knmaxtrk;k--;) {
652 dig.SetTrack(k,-3);
653 dig.SetHit(k,-1);
654 }
02d6eccc 655 aliITS->AddSimDigit(AliITSUGeomTGeo::kDetTypePix,&dig);
4fa9d550 656 if (GetDebug(2)) AliInfo(Form("Add noisy pixel %d(%d/%d) Noise=%d",ordV[ordI[j]],iz,ix,dig.GetSignalPix()));
451f5018 657 }
658 return ncand;
659}
660
661//______________________________________________________________________
662void AliITSUSimulationPix::SetCoupling(AliITSUSDigit* old, Int_t ntrack, Int_t idhit)
663{
664 // Take into account the coupling between adiacent pixels.
665 // The parameters probcol and probrow are the probability of the
666 // signal in one pixel shared in the two adjacent pixels along
667 // the column and row direction, respectively.
668 // Note pList is goten via GetMap() and module is not need any more.
669 // Otherwise it is identical to that coded by Tiziano Virgili (BSN).
670 //Begin_Html
671 /*
672 <img src="picts/ITS/barimodel_3.gif">
673 </pre>
674 <br clear=left>
675 <font size=+2 color=red>
676 <a href="mailto:tiziano.virgili@cern.ch"></a>.
677 </font>
678 <pre>
679 */
680 //End_Html
681 // Inputs:
682 // old existing AliITSUSDigit
683 // Int_t ntrack track incex number
684 // Int_t idhit hit index number
685 UInt_t col,row;
686 Double_t pulse1,pulse2;
687 Double_t couplR=0.0,couplC=0.0;
688 Double_t xr=0.;
689 //
690 fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
4fa9d550 691 fSimuParam->GetPixCouplingParam(couplC,couplR);
692 if (GetDebug(2)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d) couplC=%e couplR=%e",
693 col,row,ntrack,idhit,couplC,couplR));
451f5018 694 pulse2 = pulse1 = old->GetSignal();
4fa9d550 695 if (pulse1<fSimuParam->GetPixMinElToAdd()) return; // too small signal
451f5018 696 for (Int_t isign=-1;isign<=1;isign+=2) {
697 //
698 // loop in col direction
699 int j1 = int(col) + isign;
700 xr = gRandom->Rndm();
701 if ( !((j1<0) || (j1>fSeg->Npz()-1) || (xr>couplC)) ) UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
702 //
703 // loop in row direction
704 int j2 = int(row) + isign;
705 xr = gRandom->Rndm();
706 if ( !((j2<0) || (j2>fSeg->Npx()-1) || (xr>couplR)) ) UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
707 }
708 //
709}
710
711//______________________________________________________________________
712void AliITSUSimulationPix::SetCouplingOld(AliITSUSDigit* old, Int_t ntrack,Int_t idhit)
713{
714 // Take into account the coupling between adiacent pixels.
715 // The parameters probcol and probrow are the fractions of the
716 // signal in one pixel shared in the two adjacent pixels along
717 // the column and row direction, respectively.
718 //Begin_Html
719 /*
720 <img src="picts/ITS/barimodel_3.gif">
721 </pre>
722 <br clear=left>
723 <font size=+2 color=red>
724 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
725 </font>
726 <pre>
727 */
728 //End_Html
729 // Inputs:
730 // old existing AliITSUSDigit
731 // ntrack track incex number
732 // idhit hit index number
733 // module module number
734 //
735 UInt_t col,row;
736 Double_t pulse1,pulse2;
737 Double_t couplR=0.0,couplC=0.0;
738 //
739 fSensMap->GetMapIndex(old->GetUniqueID(),col,row);
4fa9d550 740 fSimuParam->GetPixCouplingParam(couplC,couplR);
741 if (GetDebug(3)) AliInfo(Form("(col=%d,row=%d,ntrack=%d,idhit=%d) couplC=%e couplR=%e",
742 col,row,ntrack,idhit,couplC,couplR));
743 //
744 if (old->GetSignal()<fSimuParam->GetPixMinElToAdd()) return; // too small signal
745 for (Int_t isign=-1;isign<=1;isign+=2) {// loop in col direction
746 pulse2 = pulse1 = old->GetSignal();
747 //
748 int j1 = int(col)+isign;
749 pulse1 *= couplC;
750 if ((j1<0)||(j1>fSeg->Npz()-1)||(pulse1<fSimuParam->GetPixThreshold(fModule))) pulse1 = old->GetSignal();
751 else UpdateMapSignal(UInt_t(j1),row,ntrack,idhit,pulse1);
752
753 // loop in row direction
754 int j2 = int(row) + isign;
755 pulse2 *= couplR;
756 if ((j2<0)||(j2>(fSeg->Npx()-1))||(pulse2<fSimuParam->GetPixThreshold(fModule))) pulse2 = old->GetSignal();
757 else UpdateMapSignal(col,UInt_t(j2),ntrack,idhit,pulse2);
758 } // for isign
451f5018 759}
760
761//______________________________________________________________________
762void AliITSUSimulationPix::GenerateStrobePhase()
763{
764 // Generate randomly the strobe
765 // phase w.r.t to the LHC clock
766 // Done once per event
767 const Double_t kBunchLenght = 25e-9; // LHC clock
768 fStrobePhase = ((Double_t)gRandom->Integer(fStrobeLenght))*kBunchLenght-
769 (Double_t)fStrobeLenght*kBunchLenght+
770 kBunchLenght/2;
771}
772