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