]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/AliITSUSimuParam.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
improved feeddown paramaterisation
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUSimuParam.cxx
CommitLineData
451f5018 1/**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id: AliITSUSimuParam.cxx 48165 2011-03-07 17:48:57Z masera $ */
17
18///////////////////////////////////////////////////////////////////
19// //
20// Implementation of the class to store the parameters used in //
9912c845 21// the simulation of ITS upgrade detectors //
22// //
29ad4146 23// On the fRespFunParam array content: it holds the //
24// AliITSUParamList type objects with data for response //
25// simulation of type of //
9912c845 26// detector (e.g. Class/Segmentation) used in the Config.C //
27// The convention is: //
29ad4146 28// 1) AliITSUParamList::GetUniqueID() defines detectorID //
29// (see header of the AliITSUGeomTGeo.h) for which these //
30// response data is defined //
9912c845 31// //
29ad4146 32// 2) AliITSUParamList::GetID() defines the charge spread //
33// function served by these data, for instance in case of //
34// Pixels these are the functions aliased to enums //
35// kSpreadFunGauss2D... in AliITSUSimulationPix.h //
9912c845 36// //
37// 3) Each detector class is free to interpred the content of //
29ad4146 38// AliITSUParamList. AliITSUSimulationPix, for instance requests //
39// that first AliITSUSimulationPix::kParamStart are reserved for //
40// some standard properties (like number of neighbours around the//
9912c845 41// pixel with the charge is injected to consider for the charge //
42// spread (the values may be different for different finctions) //
43// //
44// //
45// //
46// //
47// //
48// //
49// //
50// //
451f5018 51// //
52///////////////////////////////////////////////////////////////////
53#include "AliITSUSimuParam.h"
54#include "AliLog.h"
29ad4146 55#include "AliITSUParamList.h"
c92b1537 56
a11ef2e4 57using namespace TMath;
451f5018 58
59
4fa9d550 60const Double_t AliITSUSimuParam::fgkPixBiasVoltageDefault = 18.182;
29ad4146 61const Double_t AliITSUSimuParam::fgkPixThreshDefault = 20.;
62const Double_t AliITSUSimuParam::fgkPixThrSigmaDefault = 5.;
63const Double_t AliITSUSimuParam::fgkPixMinElToAddDefault = 1.;
120e5202 64const UInt_t AliITSUSimuParam::fgkPixCouplingOptDefault = AliITSUSimuParam::kNoCouplingPix;
4fa9d550 65const Double_t AliITSUSimuParam::fgkPixCouplColDefault = 0.;
66const Double_t AliITSUSimuParam::fgkPixCouplRowDefault = 0.055;
67const Double_t AliITSUSimuParam::fgkPixEccDiffDefault = 0.85;
68const Double_t AliITSUSimuParam::fgkPixLorentzHoleWeightDefault = 1.0;
451f5018 69const Double_t AliITSUSimuParam::fgkGeVtoChargeDefault = 3.6e-9;
70const Double_t AliITSUSimuParam::fgkDOverVDefault = 0.000375;
71const Double_t AliITSUSimuParam::fgkTDefault = 300;
29ad4146 72const Double_t AliITSUSimuParam::fgkPixFakeRateDefault = 1e-4;
73const Bool_t AliITSUSimuParam::fgkPixNoiseInAllMod = kFALSE;
451f5018 74
75const Double_t AliITSUSimuParam::fgkNsigmasDefault = 3.;
76const Int_t AliITSUSimuParam::fgkNcompsDefault = 121;
77
78ClassImp(AliITSUSimuParam)
79
80//______________________________________________________________________
81AliITSUSimuParam::AliITSUSimuParam()
82: fGeVcharge(fgkGeVtoChargeDefault)
83 ,fDOverV(fgkDOverVDefault)
84 ,fT(fgkTDefault)
85 //
4fa9d550 86 ,fNPix(0)
120e5202 87 ,fPixCouplOpt(kNoCouplingPix)
4fa9d550 88 ,fPixCouplCol(fgkPixCouplColDefault)
89 ,fPixCouplRow(fgkPixCouplRowDefault)
4fa9d550 90 ,fPixLorentzDrift(kTRUE)
91 ,fPixLorentzHoleWeight(fgkPixLorentzHoleWeightDefault)
92 ,fPixAddNoisyFlag(kFALSE)
93 ,fPixRemoveDeadFlag(kFALSE)
451f5018 94 //
4fa9d550 95 ,fPixThreshDef(fgkPixThreshDefault)
96 ,fPixThrSigmaDef(fgkPixThrSigmaDefault)
97 ,fPixBiasVoltageDef(fgkPixBiasVoltageDefault)
98 ,fPixNoiseDef(0)
99 ,fPixBaselineDef(0)
100 ,fPixMinElToAddDef(fgkPixMinElToAddDefault)
29ad4146 101 ,fPixFakeRateDef(fgkPixFakeRateDefault)
102 ,fPixNoiseInAllMod(fgkPixNoiseInAllMod)
451f5018 103 //
4fa9d550 104 ,fPixThresh(0)
105 ,fPixThrSigma(0)
106 ,fPixBiasVoltage(0)
107 ,fPixSigma(0)
108 ,fPixNoise(0)
109 ,fPixBaseline(0)
c92b1537 110 ,fRespFunParam(0)
451f5018 111{
112 // default constructor
c92b1537 113 fRespFunParam.SetOwner(kTRUE);
451f5018 114}
115
116//______________________________________________________________________
4fa9d550 117AliITSUSimuParam::AliITSUSimuParam(UInt_t nPix)
451f5018 118 :fGeVcharge(fgkGeVtoChargeDefault)
119 ,fDOverV(fgkDOverVDefault)
120 ,fT(fgkTDefault)
121 //
4fa9d550 122 ,fNPix(nPix)
120e5202 123 ,fPixCouplOpt(kNoCouplingPix)
4fa9d550 124 ,fPixCouplCol(fgkPixCouplColDefault)
125 ,fPixCouplRow(fgkPixCouplRowDefault)
4fa9d550 126 ,fPixLorentzDrift(kTRUE)
127 ,fPixLorentzHoleWeight(fgkPixLorentzHoleWeightDefault)
128 ,fPixAddNoisyFlag(kFALSE)
129 ,fPixRemoveDeadFlag(kFALSE)
451f5018 130 //
4fa9d550 131 ,fPixThreshDef(fgkPixThreshDefault)
132 ,fPixThrSigmaDef(fgkPixThrSigmaDefault)
133 ,fPixBiasVoltageDef(fgkPixBiasVoltageDefault)
134 ,fPixNoiseDef(0)
135 ,fPixBaselineDef(0)
136 ,fPixMinElToAddDef(fgkPixMinElToAddDefault)
29ad4146 137 ,fPixFakeRateDef(fgkPixFakeRateDefault)
138 ,fPixNoiseInAllMod(fgkPixNoiseInAllMod)
451f5018 139 //
4fa9d550 140 ,fPixThresh(0)
141 ,fPixThrSigma(0)
142 ,fPixBiasVoltage(0)
143 ,fPixSigma(0)
144 ,fPixNoise(0)
145 ,fPixBaseline(0)
c92b1537 146 ,fRespFunParam(0)
451f5018 147{
148 // regular constructor
4fa9d550 149 if (fNPix>0) {
150 fPixBiasVoltage = new Double_t[fNPix];
151 fPixThresh = new Double_t[fNPix];
152 fPixThrSigma = new Double_t[fNPix];
153 fPixNoise = new Double_t[fNPix];
154 fPixBaseline = new Double_t[fNPix];
451f5018 155 }
4fa9d550 156 SetPixThreshold(fgkPixThreshDefault,fgkPixThrSigmaDefault);
157 SetPixNoise(0.,0.);
158 SetPixBiasVoltage(fgkPixBiasVoltageDefault);
c92b1537 159 fRespFunParam.SetOwner(kTRUE);
451f5018 160 //
161}
162
163//______________________________________________________________________
164AliITSUSimuParam::AliITSUSimuParam(const AliITSUSimuParam &simpar)
165 :TObject(simpar)
166 ,fGeVcharge(simpar.fGeVcharge)
167 ,fDOverV(simpar.fDOverV)
168 ,fT(simpar.fT)
169 //
4fa9d550 170 ,fNPix(simpar.fNPix)
171 ,fPixCouplOpt(simpar.fPixCouplOpt)
172 ,fPixCouplCol(simpar.fPixCouplCol)
173 ,fPixCouplRow(simpar.fPixCouplRow)
4fa9d550 174 ,fPixLorentzDrift(simpar.fPixLorentzDrift)
175 ,fPixLorentzHoleWeight(simpar.fPixLorentzHoleWeight)
176 ,fPixAddNoisyFlag(simpar.fPixAddNoisyFlag)
177 ,fPixRemoveDeadFlag(simpar.fPixRemoveDeadFlag)
451f5018 178 //
4fa9d550 179 ,fPixThreshDef(simpar.fPixThreshDef)
180 ,fPixThrSigmaDef(simpar.fPixThrSigmaDef)
181 ,fPixBiasVoltageDef(simpar.fPixBiasVoltageDef)
182 ,fPixNoiseDef(simpar.fPixNoiseDef)
183 ,fPixBaselineDef(simpar.fPixBaselineDef)
184 ,fPixMinElToAddDef(simpar.fPixMinElToAddDef)
29ad4146 185 ,fPixFakeRateDef(simpar.fPixFakeRateDef)
186 ,fPixNoiseInAllMod(simpar.fPixNoiseInAllMod)
187 //
4fa9d550 188 ,fPixThresh(0)
189 ,fPixThrSigma(0)
190 ,fPixBiasVoltage(0)
191 ,fPixSigma(0)
192 ,fPixNoise(0)
193 ,fPixBaseline(0)
c92b1537 194 ,fRespFunParam(0)
451f5018 195 //
196{
197 // copy constructor
4fa9d550 198 if (fNPix) {
199 fPixBiasVoltage = new Double_t[fNPix];
200 fPixThresh = new Double_t[fNPix];
201 fPixThrSigma = new Double_t[fNPix];
202 fPixNoise = new Double_t[fNPix];
203 fPixBaseline = new Double_t[fNPix];
451f5018 204 }
4fa9d550 205 for (Int_t i=fNPix;i--;) {
206 fPixBiasVoltage[i] = simpar.fPixBiasVoltage[i];
207 fPixThresh[i] = simpar.fPixThresh[i];
208 fPixThrSigma[i] = simpar.fPixThrSigma[i];
209 fPixNoise[i] = simpar.fPixNoise[i];
210 fPixBaseline[i] = simpar.fPixBaseline[i];
451f5018 211 }
212 //
c92b1537 213 for (int i=0;i<simpar.fRespFunParam.GetEntriesFast();i++) {
29ad4146 214 AliITSUParamList* pr = (AliITSUParamList*)simpar.fRespFunParam[0];
215 if (pr) fRespFunParam.AddLast(new AliITSUParamList(*pr));
c92b1537 216 }
217 fRespFunParam.SetOwner(kTRUE);
451f5018 218}
219
220//______________________________________________________________________
221AliITSUSimuParam& AliITSUSimuParam::operator=(const AliITSUSimuParam& source)
222{
223 // Assignment operator.
224 if (this==&source) return *this;
225 this->~AliITSUSimuParam();
226 new(this) AliITSUSimuParam(source);
227 return *this;
228 //
229}
230
231//______________________________________________________________________
232AliITSUSimuParam::~AliITSUSimuParam()
233{
234 // destructor
4fa9d550 235 delete[] fPixBiasVoltage;
236 delete[] fPixThresh;
237 delete[] fPixThrSigma;
238 delete[] fPixNoise;
239 delete[] fPixBaseline;
451f5018 240}
241
242//________________________________________________________________________
243void AliITSUSimuParam::Print(Option_t *) const
244{
245 // Dump all parameters
246 Dump();
29ad4146 247 //
248 int nresp = fRespFunParam.GetEntriesFast();
249 if (nresp) {
250 printf("Individual sensor responses\n");
251 for (int i=0;i<nresp;i++) {
252 AliITSUParamList* lst = (AliITSUParamList*)fRespFunParam[i];
253 if (!lst) continue;
254 lst->Print();
255 }
256 }
451f5018 257}
258
259//_______________________________________________________________________
4fa9d550 260Double_t AliITSUSimuParam::ApplyPixBaselineAndNoise(UInt_t mod) const
451f5018 261{
262 // generate random noise
263 double base,noise;
4fa9d550 264 if (mod>=fNPix) {
265 if (fNPix>0) {AliFatal(Form("Wrong module %d, NPidUpg=%d",mod,fNPix));}
266 base = fPixBaselineDef;
267 noise = fPixNoiseDef;
451f5018 268 }
269 else {
4fa9d550 270 base = fPixBaseline[mod];
271 noise = fPixNoise[mod];
451f5018 272 }
273 return base+noise*gRandom->Gaus();
274}
275
276//_______________________________________________________________________
277Double_t AliITSUSimuParam::CalcProbNoiseOverThreshold(UInt_t mod) const
278{
279 // calculate probability of noise exceeding the threshold
280 double base,noise,thresh;
4fa9d550 281 if (mod>=fNPix) {
282 if (fNPix>0) {AliFatal(Form("Wrong module %d, NPidUpg=%d",mod,fNPix));}
283 base = fPixBaselineDef;
284 noise = fPixNoiseDef;
285 thresh = fPixThreshDef;
451f5018 286 }
287 else {
4fa9d550 288 base = fPixBaseline[mod];
289 noise = fPixNoise[mod];
290 thresh = fPixThresh[mod];
451f5018 291 }
292 if (noise<1e-12) {
293 if (base>thresh) return 1;
294 else return 0;
295 }
296 return CalcProbNoiseOverThreshold(base, noise, thresh);
297}
298
299//_______________________________________________________________________
4fa9d550 300void AliITSUSimuParam::SetPixThreshold(Double_t thresh, Double_t sigma, int mod)
451f5018 301{
302 // set threshold params
303 if (mod<0) {
4fa9d550 304 fPixThreshDef = thresh;
305 fPixThrSigmaDef = sigma;
306 for (int i=fNPix;i--;) {
307 fPixThresh[i] = thresh;
308 fPixThrSigma[i] = sigma;
451f5018 309 }
310 }
4fa9d550 311 else if (mod>=(int)fNPix) {
312 if (fNPix>0) {AliFatal(Form("Wrong module %d, NPidUpg=%d",mod,fNPix));}
313 fPixThreshDef = thresh;
314 fPixThrSigmaDef = sigma;
451f5018 315 }
316 else {
4fa9d550 317 fPixThresh[mod] = thresh;
318 fPixThrSigma[mod] = sigma;
451f5018 319 }
320 //
321}
322
323//_______________________________________________________________________
4fa9d550 324Double_t AliITSUSimuParam::GetPixThreshold(UInt_t mod) const
451f5018 325{
326 // obtain threshold
4fa9d550 327 if (mod>=fNPix) {
328 if (fNPix>0) {AliFatal(Form("Wrong module %d, NPidUpg=%d",mod,fNPix));}
329 return fPixThreshDef;
451f5018 330 }
4fa9d550 331 else return fPixThresh[mod];
451f5018 332}
333
334//_______________________________________________________________________
4fa9d550 335void AliITSUSimuParam::GetPixThreshold(UInt_t mod, Double_t &thresh, Double_t &sigma) const
451f5018 336{
337 // obtain thresholds
4fa9d550 338 if (mod>=fNPix) {
339 if (fNPix>0) {AliFatal(Form("Wrong module %d, NPidUpg=%d",mod,fNPix));}
340 thresh = fPixThreshDef;
341 sigma = fPixThrSigmaDef;
451f5018 342 }
343 else {
4fa9d550 344 thresh = fPixThresh[mod];
345 sigma = fPixThrSigma[mod];
451f5018 346 }
347}
348
349//_______________________________________________________________________
4fa9d550 350void AliITSUSimuParam::SetPixBiasVoltage(Double_t val, int mod)
451f5018 351{
352 // set threshold params
353 if (mod<0) {
4fa9d550 354 fPixBiasVoltageDef = val;
355 for (int i=fNPix;i--;) fPixBiasVoltage[i] = val;
451f5018 356 }
4fa9d550 357 else if (mod>=(int)fNPix) {
358 if (fNPix>0) {AliFatal(Form("Wrong module %d, NPidUpg=%d",mod,fNPix));}
359 fPixBiasVoltageDef = val;
451f5018 360 }
4fa9d550 361 else fPixBiasVoltage[mod] = val;
451f5018 362 //
363}
364
365//_______________________________________________________________________
4fa9d550 366Double_t AliITSUSimuParam::GetPixBiasVoltage(UInt_t mod) const
451f5018 367{
368 // obtain threshold
4fa9d550 369 if (mod>=fNPix) {
370 if (fNPix>0) {AliFatal(Form("Wrong module %d, NPidUpg=%d",mod,fNPix));}
371 return fPixBiasVoltageDef;
451f5018 372 }
4fa9d550 373 else return fPixBiasVoltage[mod];
451f5018 374}
375
376//_______________________________________________________________________
4fa9d550 377void AliITSUSimuParam::SetPixNoise(Double_t noise, Double_t baseline, int mod)
451f5018 378{
379 // set noise params
380 if (mod<0) {
4fa9d550 381 fPixNoiseDef = noise;
382 fPixBaselineDef = baseline;
383 for (int i=fNPix;i--;) {
384 fPixNoise[i] = noise;
385 fPixBaseline[i] = baseline;
451f5018 386 }
387 }
4fa9d550 388 else if (mod>=(int)fNPix) {
389 if (fNPix>0) {AliFatal(Form("Wrong module %d, NPidUpg=%d",mod,fNPix));}
390 fPixNoiseDef = noise;
391 fPixBaselineDef = baseline;
451f5018 392 }
393 else {
4fa9d550 394 fPixNoise[mod] = noise;
395 fPixBaseline[mod] = baseline;
451f5018 396 }
397 //
398}
399
400//_______________________________________________________________________
4fa9d550 401void AliITSUSimuParam::GetPixNoise(UInt_t mod, Double_t &noise, Double_t &baseline) const
451f5018 402{
403 // obtain noise
4fa9d550 404 if (mod>=fNPix) {
405 if (fNPix>0) {AliFatal(Form("Wrong module %d, NPidUpg=%d",mod,fNPix));}
406 noise = fPixNoiseDef;
407 baseline = fPixBaselineDef;
451f5018 408 }
409 else {
4fa9d550 410 noise = fPixNoise[mod];
411 baseline = fPixBaseline[mod];
451f5018 412 }
413}
414
415//_______________________________________________________________________
4fa9d550 416void AliITSUSimuParam::SetPixCouplingOption(UInt_t opt)
451f5018 417{
418 // set coupling option
4fa9d550 419 if (opt>=kMaxCouplingOptPix) AliFatal(Form("Coupling option %d should be less than %d",opt,kMaxCouplingOptPix));
420 fPixCouplOpt = opt;
451f5018 421}
422
423
424//______________________________________________________________________
425Double_t AliITSUSimuParam::LorentzAngleHole(Double_t B) const
426{
427 // Computes the Lorentz angle for electrons in Si
428 // Input: magnetic Field in KGauss
429 // Output: Lorentz angle in radians (positive if Bz is positive)
430 // Main Reference: NIM A 497 (2003) 389–396.
431 // "An algorithm for calculating the Lorentz angle in silicon detectors", V. Bartsch et al.
432 //
433 const Double_t krH=0.70; // Hall scattering factor for Hole
434 const Double_t kT0 = 300.; // reference Temperature (degree K).
435 const Double_t kmulow0 = 470.5; // cm^2/Volt-sec
436 const Double_t keT0 = -2.5; // Power of Temp.
437 const Double_t beta0 = 1.213; // beta coeff. at T0=300K
438 const Double_t keT1 = 0.17; // Power of Temp. for beta
439 const Double_t kvsat0 = 8.37E+06; // saturated velocity at T0=300K (cm/sec)
440 const Double_t keT2 = 0.52; // Power of Temp. for vsat
441 Double_t tT = fT;
442 Double_t eE= 1./fDOverV;
a11ef2e4 443 Double_t muLow=kmulow0*Power(tT/kT0,keT0);
444 Double_t beta=beta0*Power(tT/kT0,keT1);
445 Double_t vsat=kvsat0*Power(tT/kT0,keT2);
446 Double_t mu=muLow/Power(1+Power(muLow*eE/vsat,beta),1/beta);
447 Double_t angle=ATan(krH*mu*B*1.E-05); // Conversion Factor
451f5018 448 return angle;
449}
450
451//______________________________________________________________________
452Double_t AliITSUSimuParam::LorentzAngleElectron(Double_t B) const
453{
454 // Computes the Lorentz angle for electrons in Si
455 // Input: magnetic Field in KGauss
456 // Output: Lorentz angle in radians (positive if Bz is positive)
457 // Main Reference: NIM A 497 (2003) 389–396.
458 // "An algorithm for calculating the Lorentz angle in silicon detectors", V. Bartsch et al.
459 //
460 const Double_t krH=1.15; // Hall scattering factor for Electron
461 const Double_t kT0 = 300.; // reference Temperature (degree K).
462 const Double_t kmulow0 = 1417.0; // cm^2/Volt-sec
463 const Double_t keT0 = -2.2; // Power of Temp.
464 const Double_t beta0 = 1.109; // beta coeff. at T0=300K
465 const Double_t keT1 = 0.66; // Power of Temp. for beta
466 const Double_t kvsat0 = 1.07E+07; // saturated velocity at T0=300K (cm/sec)
467 const Double_t keT2 = 0.87; // Power of Temp. for vsat
468 Double_t tT = fT;
469 Double_t eE= 1./fDOverV;
a11ef2e4 470 Double_t muLow=kmulow0*Power(tT/kT0,keT0);
471 Double_t beta=beta0*Power(tT/kT0,keT1);
472 Double_t vsat=kvsat0*Power(tT/kT0,keT2);
473 Double_t mu=muLow/Power(1+Power(muLow*eE/vsat,beta),1/beta);
474 Double_t angle=ATan(krH*mu*B*1.E-05);
451f5018 475 return angle;
476}
477
c92b1537 478//___________________________________________________________
29ad4146 479const AliITSUParamList* AliITSUSimuParam::FindRespFunParams(Int_t detId) const
451f5018 480{
c92b1537 481 // find parameters list for detID
482 for (int i=fRespFunParam.GetEntriesFast();i--;) {
29ad4146 483 const AliITSUParamList* pr = GetRespFunParams(i);
c92b1537 484 if (int(pr->GetUniqueID())==detId) return pr;
485 }
486 return 0;
451f5018 487}
488
c92b1537 489//___________________________________________________________
29ad4146 490void AliITSUSimuParam::AddRespFunParam(AliITSUParamList* pr)
451f5018 491{
c92b1537 492 // add spread parameterization data
493 fRespFunParam.AddLast(pr);
451f5018 494}
29ad4146 495//______________________________________________________________________________________
496Double_t AliITSUSimuParam::CalcProbNoiseOverThreshold(double mean, double sigma, double thresh)
497{
498 // calculate probability of noise exceeding the threshold
499 //if (mean+6*sigma<thresh) return 0;
500 //if (mean-6*sigma>thresh) return 1.;
501 //const double ksqrt2 = 1.41421356237309515e+00;
502 //return 0.5*AliMathBase::ErfcFast( (thresh-mean)/(sigma*ksqrt2));
503
504 // The MIMOSA sensors show Landau noise... RTF
505 // Calculate probability of Landau noise exceeding the threshold
506
507
508 //Double_t prob = this->LandauCdf(thresh,sigma,mean);// myLcdf->Eval(thresh);
509 //
510 //printf("====> After LandauCdf(thresh,sigma,mean)");
511
512 //return prob;
513 // (double x, double xi, double x0) {
514
515 // Implementation is taken from ROOT
516// Printf("AliITSUSimuParam::LandauCdf");
517 //arguments threshold, sigma, mean
518
519 double x = thresh;
520 double xi = sigma;
521 double x0 = mean;
522
523
524 static double p1[5] = {0.2514091491e+0,-0.6250580444e-1, 0.1458381230e-1, -0.2108817737e-2, 0.7411247290e-3};
525 static double q1[5] = {1.0 ,-0.5571175625e-2, 0.6225310236e-1, -0.3137378427e-2, 0.1931496439e-2};
526
527 static double p2[4] = {0.2868328584e+0, 0.3564363231e+0, 0.1523518695e+0, 0.2251304883e-1};
528 static double q2[4] = {1.0 , 0.6191136137e+0, 0.1720721448e+0, 0.2278594771e-1};
529
530 static double p3[4] = {0.2868329066e+0, 0.3003828436e+0, 0.9950951941e-1, 0.8733827185e-2};
531 static double q3[4] = {1.0 , 0.4237190502e+0, 0.1095631512e+0, 0.8693851567e-2};
532
533 static double p4[4] = {0.1000351630e+1, 0.4503592498e+1, 0.1085883880e+2, 0.7536052269e+1};
534 static double q4[4] = {1.0 , 0.5539969678e+1, 0.1933581111e+2, 0.2721321508e+2};
535
536 static double p5[4] = {0.1000006517e+1, 0.4909414111e+2, 0.8505544753e+2, 0.1532153455e+3};
537 static double q5[4] = {1.0 , 0.5009928881e+2, 0.1399819104e+3, 0.4200002909e+3};
538
539 static double p6[4] = {0.1000000983e+1, 0.1329868456e+3, 0.9162149244e+3, -0.9605054274e+3};
540 static double q6[4] = {1.0 , 0.1339887843e+3, 0.1055990413e+4, 0.5532224619e+3};
541
542 static double a1[4] = {0, -0.4583333333e+0, 0.6675347222e+0,-0.1641741416e+1};
543
544 static double a2[4] = {0, 1.0 ,-0.4227843351e+0,-0.2043403138e+1};
545
546 double v = (x - x0)/xi;
547 double u;
548 double lan;
549
550
551
552 if (v < -5.5) {
553 u = std::exp(v+1);
554 lan = 0.3989422803*std::exp(-1./u)*std::sqrt(u)*(1+(a1[1]+(a1[2]+a1[3]*u)*u)*u);
555 }
556 else if (v < -1 ) {
557 u = std::exp(-v-1);
558 lan = (std::exp(-u)/std::sqrt(u))*(p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v)*v)*v)*v)/
559 (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v)*v)*v)*v);
560 }
561 else if (v < 1)
562 lan = (p2[0]+(p2[1]+(p2[2]+p2[3]*v)*v)*v)/(q2[0]+(q2[1]+(q2[2]+q2[3]*v)*v)*v);
563 else if (v < 4)
564 lan = (p3[0]+(p3[1]+(p3[2]+p3[3]*v)*v)*v)/(q3[0]+(q3[1]+(q3[2]+q3[3]*v)*v)*v);
565 else if (v < 12) {
566 u = 1./v;
567 lan = (p4[0]+(p4[1]+(p4[2]+p4[3]*u)*u)*u)/(q4[0]+(q4[1]+(q4[2]+q4[3]*u)*u)*u);
568 }
569 else if (v < 50) {
570 u = 1./v;
571 lan = (p5[0]+(p5[1]+(p5[2]+p5[3]*u)*u)*u)/(q5[0]+(q5[1]+(q5[2]+q5[3]*u)*u)*u);
572 }
573 else if (v < 300) {
574 u = 1./v;
575 lan = (p6[0]+(p6[1]+(p6[2]+p6[3]*u)*u)*u)/(q6[0]+(q6[1]+(q6[2]+q6[3]*u)*u)*u);
576 }
577 else {
578 u = 1./(v-v*std::log(v)/(v+1));
579 lan = 1-(a2[1]+(a2[2]+a2[3]*u)*u)*u;
580 }
581
582
583 return lan;
584
585
586
587
588}
589
590//_______________________________________________________________________
591Double_t AliITSUSimuParam::GenerateNoiseQFunction(double prob, double mean, double sigma)
592{
593 // generate random noise exceeding threshold probability prob, i.e. find a random point in the right
594 // tail of the gaussian(base,noise), provided that the tail integral = prob
595 //const double ksqrt2 = 1.41421356237309515e+00;
596 //return mean+sigma*ksqrt2*TMath::ErfcInverse(2*prob*(1.-gRandom->Rndm()));
597
598 double z = gRandom->Uniform(prob,(1.0-1e-9));
599 double xi = sigma;
600
601 static double f[982] = {
602 0 , 0 , 0 ,0 ,0 ,-2.244733,
603 -2.204365,-2.168163,-2.135219,-2.104898,-2.076740,-2.050397,
604 -2.025605,-2.002150,-1.979866,-1.958612,-1.938275,-1.918760,
605 -1.899984,-1.881879,-1.864385,-1.847451,-1.831030,-1.815083,
606 -1.799574,-1.784473,-1.769751,-1.755383,-1.741346,-1.727620,
607 -1.714187,-1.701029,-1.688130,-1.675477,-1.663057,-1.650858,
608 -1.638868,-1.627078,-1.615477,-1.604058,-1.592811,-1.581729,
609 -1.570806,-1.560034,-1.549407,-1.538919,-1.528565,-1.518339,
610 -1.508237,-1.498254,-1.488386,-1.478628,-1.468976,-1.459428,
611 -1.449979,-1.440626,-1.431365,-1.422195,-1.413111,-1.404112,
612 -1.395194,-1.386356,-1.377594,-1.368906,-1.360291,-1.351746,
613 -1.343269,-1.334859,-1.326512,-1.318229,-1.310006,-1.301843,
614 -1.293737,-1.285688,-1.277693,-1.269752,-1.261863,-1.254024,
615 -1.246235,-1.238494,-1.230800,-1.223153,-1.215550,-1.207990,
616 -1.200474,-1.192999,-1.185566,-1.178172,-1.170817,-1.163500,
617 -1.156220,-1.148977,-1.141770,-1.134598,-1.127459,-1.120354,
618 -1.113282,-1.106242,-1.099233,-1.092255,
619 -1.085306,-1.078388,-1.071498,-1.064636,-1.057802,-1.050996,
620 -1.044215,-1.037461,-1.030733,-1.024029,-1.017350,-1.010695,
621 -1.004064, -.997456, -.990871, -.984308, -.977767, -.971247,
622 -.964749, -.958271, -.951813, -.945375, -.938957, -.932558,
623 -.926178, -.919816, -.913472, -.907146, -.900838, -.894547,
624 -.888272, -.882014, -.875773, -.869547, -.863337, -.857142,
625 -.850963, -.844798, -.838648, -.832512, -.826390, -.820282,
626 -.814187, -.808106, -.802038, -.795982, -.789940, -.783909,
627 -.777891, -.771884, -.765889, -.759906, -.753934, -.747973,
628 -.742023, -.736084, -.730155, -.724237, -.718328, -.712429,
629 -.706541, -.700661, -.694791, -.688931, -.683079, -.677236,
630 -.671402, -.665576, -.659759, -.653950, -.648149, -.642356,
631 -.636570, -.630793, -.625022, -.619259, -.613503, -.607754,
632 -.602012, -.596276, -.590548, -.584825, -.579109, -.573399,
633 -.567695, -.561997, -.556305, -.550618, -.544937, -.539262,
634 -.533592, -.527926, -.522266, -.516611, -.510961, -.505315,
635 -.499674, -.494037, -.488405, -.482777,
636 -.477153, -.471533, -.465917, -.460305, -.454697, -.449092,
637 -.443491, -.437893, -.432299, -.426707, -.421119, -.415534,
638 -.409951, -.404372, -.398795, -.393221, -.387649, -.382080,
639 -.376513, -.370949, -.365387, -.359826, -.354268, -.348712,
640 -.343157, -.337604, -.332053, -.326503, -.320955, -.315408,
641 -.309863, -.304318, -.298775, -.293233, -.287692, -.282152,
642 -.276613, -.271074, -.265536, -.259999, -.254462, -.248926,
643 -.243389, -.237854, -.232318, -.226783, -.221247, -.215712,
644 -.210176, -.204641, -.199105, -.193568, -.188032, -.182495,
645 -.176957, -.171419, -.165880, -.160341, -.154800, -.149259,
646 -.143717, -.138173, -.132629, -.127083, -.121537, -.115989,
647 -.110439, -.104889, -.099336, -.093782, -.088227, -.082670,
648 -.077111, -.071550, -.065987, -.060423, -.054856, -.049288,
649 -.043717, -.038144, -.032569, -.026991, -.021411, -.015828,
650 -.010243, -.004656, .000934, .006527, .012123, .017722,
651 .023323, .028928, .034535, .040146, .045759, .051376,
652 .056997, .062620, .068247, .073877,
653 .079511, .085149, .090790, .096435, .102083, .107736,
654 .113392, .119052, .124716, .130385, .136057, .141734,
655 .147414, .153100, .158789, .164483, .170181, .175884,
656 .181592, .187304, .193021, .198743, .204469, .210201,
657 .215937, .221678, .227425, .233177, .238933, .244696,
658 .250463, .256236, .262014, .267798, .273587, .279382,
659 .285183, .290989, .296801, .302619, .308443, .314273,
660 .320109, .325951, .331799, .337654, .343515, .349382,
661 .355255, .361135, .367022, .372915, .378815, .384721,
662 .390634, .396554, .402481, .408415, .414356, .420304,
663 .426260, .432222, .438192, .444169, .450153, .456145,
664 .462144, .468151, .474166, .480188, .486218, .492256,
665 .498302, .504356, .510418, .516488, .522566, .528653,
666 .534747, .540850, .546962, .553082, .559210, .565347,
667 .571493, .577648, .583811, .589983, .596164, .602355,
668 .608554, .614762, .620980, .627207, .633444, .639689,
669 .645945, .652210, .658484, .664768,
670 .671062, .677366, .683680, .690004, .696338, .702682,
671 .709036, .715400, .721775, .728160, .734556, .740963,
672 .747379, .753807, .760246, .766695, .773155, .779627,
673 .786109, .792603, .799107, .805624, .812151, .818690,
674 .825241, .831803, .838377, .844962, .851560, .858170,
675 .864791, .871425, .878071, .884729, .891399, .898082,
676 .904778, .911486, .918206, .924940, .931686, .938446,
677 .945218, .952003, .958802, .965614, .972439, .979278,
678 .986130, .992996, .999875, 1.006769, 1.013676, 1.020597,
679 1.027533, 1.034482, 1.041446, 1.048424, 1.055417, 1.062424,
680 1.069446, 1.076482, 1.083534, 1.090600, 1.097681, 1.104778,
681 1.111889, 1.119016, 1.126159, 1.133316, 1.140490, 1.147679,
682 1.154884, 1.162105, 1.169342, 1.176595, 1.183864, 1.191149,
683 1.198451, 1.205770, 1.213105, 1.220457, 1.227826, 1.235211,
684 1.242614, 1.250034, 1.257471, 1.264926, 1.272398, 1.279888,
685 1.287395, 1.294921, 1.302464, 1.310026, 1.317605, 1.325203,
686 1.332819, 1.340454, 1.348108, 1.355780,
687 1.363472, 1.371182, 1.378912, 1.386660, 1.394429, 1.402216,
688 1.410024, 1.417851, 1.425698, 1.433565, 1.441453, 1.449360,
689 1.457288, 1.465237, 1.473206, 1.481196, 1.489208, 1.497240,
690 1.505293, 1.513368, 1.521465, 1.529583, 1.537723, 1.545885,
691 1.554068, 1.562275, 1.570503, 1.578754, 1.587028, 1.595325,
692 1.603644, 1.611987, 1.620353, 1.628743, 1.637156, 1.645593,
693 1.654053, 1.662538, 1.671047, 1.679581, 1.688139, 1.696721,
694 1.705329, 1.713961, 1.722619, 1.731303, 1.740011, 1.748746,
695 1.757506, 1.766293, 1.775106, 1.783945, 1.792810, 1.801703,
696 1.810623, 1.819569, 1.828543, 1.837545, 1.846574, 1.855631,
697 1.864717, 1.873830, 1.882972, 1.892143, 1.901343, 1.910572,
698 1.919830, 1.929117, 1.938434, 1.947781, 1.957158, 1.966566,
699 1.976004, 1.985473, 1.994972, 2.004503, 2.014065, 2.023659,
700 2.033285, 2.042943, 2.052633, 2.062355, 2.072110, 2.081899,
701 2.091720, 2.101575, 2.111464, 2.121386, 2.131343, 2.141334,
702 2.151360, 2.161421, 2.171517, 2.181648, 2.191815, 2.202018,
703 2.212257, 2.222533, 2.232845, 2.243195,
704 2.253582, 2.264006, 2.274468, 2.284968, 2.295507, 2.306084,
705 2.316701, 2.327356, 2.338051, 2.348786, 2.359562, 2.370377,
706 2.381234, 2.392131, 2.403070, 2.414051, 2.425073, 2.436138,
707 2.447246, 2.458397, 2.469591, 2.480828, 2.492110, 2.503436,
708 2.514807, 2.526222, 2.537684, 2.549190, 2.560743, 2.572343,
709 2.583989, 2.595682, 2.607423, 2.619212, 2.631050, 2.642936,
710 2.654871, 2.666855, 2.678890, 2.690975, 2.703110, 2.715297,
711 2.727535, 2.739825, 2.752168, 2.764563, 2.777012, 2.789514,
712 2.802070, 2.814681, 2.827347, 2.840069, 2.852846, 2.865680,
713 2.878570, 2.891518, 2.904524, 2.917588, 2.930712, 2.943894,
714 2.957136, 2.970439, 2.983802, 2.997227, 3.010714, 3.024263,
715 3.037875, 3.051551, 3.065290, 3.079095, 3.092965, 3.106900,
716 3.120902, 3.134971, 3.149107, 3.163312, 3.177585, 3.191928,
717 3.206340, 3.220824, 3.235378, 3.250005, 3.264704, 3.279477,
718 3.294323, 3.309244, 3.324240, 3.339312, 3.354461, 3.369687,
719 3.384992, 3.400375, 3.415838, 3.431381, 3.447005, 3.462711,
720 3.478500, 3.494372, 3.510328, 3.526370,
721 3.542497, 3.558711, 3.575012, 3.591402, 3.607881, 3.624450,
722 3.641111, 3.657863, 3.674708, 3.691646, 3.708680, 3.725809,
723 3.743034, 3.760357, 3.777779, 3.795300, 3.812921, 3.830645,
724 3.848470, 3.866400, 3.884434, 3.902574, 3.920821, 3.939176,
725 3.957640, 3.976215, 3.994901, 4.013699, 4.032612, 4.051639,
726 4.070783, 4.090045, 4.109425, 4.128925, 4.148547, 4.168292,
727 4.188160, 4.208154, 4.228275, 4.248524, 4.268903, 4.289413,
728 4.310056, 4.330832, 4.351745, 4.372794, 4.393982, 4.415310,
729 4.436781, 4.458395, 4.480154, 4.502060, 4.524114, 4.546319,
730 4.568676, 4.591187, 4.613854, 4.636678, 4.659662, 4.682807,
731 4.706116, 4.729590, 4.753231, 4.777041, 4.801024, 4.825179,
732 4.849511, 4.874020, 4.898710, 4.923582, 4.948639, 4.973883,
733 4.999316, 5.024942, 5.050761, 5.076778, 5.102993, 5.129411,
734 5.156034, 5.182864, 5.209903, 5.237156, 5.264625, 5.292312,
735 5.320220, 5.348354, 5.376714, 5.405306, 5.434131, 5.463193,
736 5.492496, 5.522042, 5.551836, 5.581880, 5.612178, 5.642734,
737 5.673552, 5.704634, 5.735986, 5.767610,
738 5.799512, 5.831694, 5.864161, 5.896918, 5.929968, 5.963316,
739 5.996967, 6.030925, 6.065194, 6.099780, 6.134687, 6.169921,
740 6.205486, 6.241387, 6.277630, 6.314220, 6.351163, 6.388465,
741 6.426130, 6.464166, 6.502578, 6.541371, 6.580553, 6.620130,
742 6.660109, 6.700495, 6.741297, 6.782520, 6.824173, 6.866262,
743 6.908795, 6.951780, 6.995225, 7.039137, 7.083525, 7.128398,
744 7.173764, 7.219632, 7.266011, 7.312910, 7.360339, 7.408308,
745 7.456827, 7.505905, 7.555554, 7.605785, 7.656608, 7.708035,
746 7.760077, 7.812747, 7.866057, 7.920019, 7.974647, 8.029953,
747 8.085952, 8.142657, 8.200083, 8.258245, 8.317158, 8.376837,
748 8.437300, 8.498562, 8.560641, 8.623554, 8.687319, 8.751955,
749 8.817481, 8.883916, 8.951282, 9.019600, 9.088889, 9.159174,
750 9.230477, 9.302822, 9.376233, 9.450735, 9.526355, 9.603118,
751 9.681054, 9.760191, 9.840558, 9.922186,10.005107,10.089353,
752 10.174959,10.261958,10.350389,10.440287,10.531693,10.624646,
753 10.719188,10.815362,10.913214,11.012789,11.114137,11.217307,
754 11.322352,11.429325,11.538283,11.649285,
755 11.762390,11.877664,11.995170,12.114979,12.237161,12.361791,
756 12.488946,12.618708,12.751161,12.886394,13.024498,13.165570,
757 13.309711,13.457026,13.607625,13.761625,13.919145,14.080314,
758 14.245263,14.414134,14.587072,14.764233,14.945778,15.131877,
759 15.322712,15.518470,15.719353,15.925570,16.137345,16.354912,
760 16.578520,16.808433,17.044929,17.288305,17.538873,17.796967,
761 18.062943,18.337176,18.620068,18.912049,19.213574,19.525133,
762 19.847249,20.180480,20.525429,20.882738,21.253102,21.637266,
763 22.036036,22.450278,22.880933,23.329017,23.795634,24.281981,
764 24.789364,25.319207,25.873062,26.452634,27.059789,27.696581,
765 28.365274,29.068370,29.808638,30.589157,31.413354,32.285060,
766 33.208568,34.188705,35.230920,36.341388,37.527131,38.796172,
767 40.157721,41.622399,43.202525,44.912465,46.769077,48.792279,
768 51.005773,53.437996,56.123356,59.103894 };
769
770 if (xi <= 0) return 0;
771 // if (z <= 0) return -std::numeric_limits<double>::infinity();
772 //if (z >= 1) return std::numeric_limits<double>::infinity();
773
774 double ranlan, u, v;
775 u = 1000*z;
776 int i = int(u);
777 u -= i;
778 if (i >= 70 && i < 800) {
779 ranlan = f[i-1] + u*(f[i] - f[i-1]);
780 } else if (i >= 7 && i <= 980) {
781 ranlan = f[i-1] + u*(f[i]-f[i-1]-0.25*(1-u)*(f[i+1]-f[i]-f[i-1]+f[i-2]));
782 } else if (i < 7) {
783 v = std::log(z);
784 u = 1/v;
785 ranlan = ((0.99858950+(3.45213058E1+1.70854528E1*u)*u)/
786 (1 +(3.41760202E1+4.01244582 *u)*u))*
787 (-std::log(-0.91893853-v)-1);
788 } else {
789 u = 1-z;
790 v = u*u;
791 if (z <= 0.999) {
792 ranlan = (1.00060006+2.63991156E2*u+4.37320068E3*v)/
793 ((1 +2.57368075E2*u+3.41448018E3*v)*u);
794 } else {
795 ranlan = (1.00001538+6.07514119E3*u+7.34266409E5*v)/
796 ((1 +6.06511919E3*u+6.94021044E5*v)*u);
797 }
798 }
799
800
801
802 Double_t lq = xi*ranlan;
803
804
805 Double_t qOverThresh = lq + mean;
806
807 return qOverThresh;
808
809
810}
811//_______________________________________________________________________
812
813
814
815