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