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