65 #if !defined(__CINT__)
68 #if !defined(__CINT__) || defined(__MAKECINT__)
75 #if !defined(__CINT__) || defined(__MAKECINT__)
76 #include "Riostream.h"
82 #include "TClassTable.h"
83 #include "StBichsel/Bichsel.h"
84 #include "StBichsel/StdEdxModel.h"
85 #include "StBichsel/StdEdxPull.h"
95 double operator() (
double *x,
double * )
const {
96 return fFunc->Derivative(*x);
101 TF1 *Derivative(TF1 *f1) {
103 TString name(
"der_");
104 name += f1->GetName();
105 TF1 * f2 =
new TF1(name,deriv, f1->GetXmin(), f1->GetXmax(), 0,
"MyDerivFunc");
109 const Int_t NMasses = 20;
110 const Double_t kAu2Gev=0.9314943228;
118 const Int_t kpionIdx = 2;
121 {
"e", 0, 1, 0, 0.51099907e-3},
122 {
"#mu", 4, 1, 1, 0.1056584},
123 {
"#pi", 3, 1, 2, 0.13956995},
124 {
"K", 2, 1, 3, 0.493677},
125 {
"p", 1, 1, 4, 0.93827231},
126 {
"d", 5, 1, 5, 1.87561339},
127 {
"t", 6, 1, 6, 2.80925},
128 {
"He3", 7, 2, 7, 2.80923},
129 {
"#alpha", 8, 2, 8, 4.00260325415*kAu2Gev},
130 {
"He6", -1, 2, 12, 6.018889124*kAu2Gev},
131 {
"He8", -1, 2, 12, 8.033921897*kAu2Gev},
132 {
"Li6", 10, 3, 10, 5.6031},
133 {
"Li7", 9, 3, 9, 6.5354},
134 {
"Be7", -1, 4, 0, 7.016003437*kAu2Gev},
135 {
"Be9", -1, 4, 0, 9.01218307*kAu2Gev},
136 {
"Be10", -1, 4, 0, 10.01353470*kAu2Gev},
137 {
"B10", -1, 5, 0, 10.01353470*kAu2Gev},
138 {
"B11", -1, 5, 0, 11.009305167*kAu2Gev},
139 {
"2#pi", 3, -1, -2, -0.13956995},
140 {
"2p", 1, -1, 0, -0.93827231}
143 const Char_t *FNames[11] = {
"Girrf",
"Sirrf",
"z",
"I70",
"I60",
"I70M",
"dNdx",
"zM",
"zN",
"70Trs",
"zTrs"};
144 const Int_t Nlog2dx = 3;
145 const Double_t log2dx[Nlog2dx] = {0,1,2};
146 static Bool_t fgRigidity = kFALSE;
148 Double_t pOverM(Double_t x) {
149 if (! fgRigidity)
return TMath::Power(10.,x);
150 return TMath::Power(10.,TMath::Abs(x) - 1.5);
153 Double_t bichselZ(Double_t *x,Double_t *par) {
154 Double_t pove = pOverM(x[0]);
156 Double_t mass = par[0];
157 if (mass < 0) {mass = - mass; scale = 2;}
158 Double_t poverm = pove/mass;
159 Double_t charge = 1.;
164 dx2 = TMath::Log2(5.);
166 return TMath::Log10(scale*charge*charge*TMath::Exp(Bichsel::Instance()->GetMostProbableZ(TMath::Log10(poverm),dx2)));
171 Double_t dEdxModelZ(Double_t *x,Double_t *par) {
172 Double_t pove = pOverM(x[0]);
174 Double_t mass = par[0];
175 if (mass < 0) {mass = - mass; scale = 2;}
176 Double_t poverm = pove/mass;
177 Double_t charge = par[1];
178 poverm *= TMath::Abs(charge);
179 return TMath::Log10(1e6*scale*StdEdxPull::EvalPred(poverm, 1, charge));
182 Double_t bichselZM(Double_t *x,Double_t *par) {
183 Double_t pove = pOverM(x[0]);
185 Double_t mass = par[0];
186 if (mass < 0) {mass = - mass; scale = 2;}
187 Double_t poverm = pove/mass;
188 Double_t charge = 1.;
192 poverm *= TMath::Abs(charge);
193 dx2 = TMath::Log2(5.);
195 return TMath::Log10(1e6*scale*StdEdxPull::EvalPred(poverm, 1, charge));
198 Double_t bichsel70(Double_t *x,Double_t *par) {
199 Double_t pove = pOverM(x[0]);
201 Double_t mass = par[0];
202 if (mass < 0) {mass = - mass; scale = 2;}
203 Double_t poverm = pove/mass;
204 Double_t charge = 1.;
209 dx2 = TMath::Log2(5.);
211 return TMath::Log10(1e6*scale*StdEdxPull::EvalPred(poverm, 0, charge));
214 Double_t bichsel70M(Double_t *x,Double_t *par) {
215 Double_t pove = pOverM(x[0]);
217 Double_t mass = par[0];
218 if (mass < 0) {mass = - mass; scale = 2;}
219 Double_t poverm = pove/mass;
220 Double_t charge = 1.;
225 dx2 = TMath::Log2(5.);
227 return TMath::Log10(1e6*scale*StdEdxPull::EvalPred(poverm, 0, charge));
231 Double_t bichsel70Trs(Double_t *x,Double_t *par) {
232 Double_t pove = pOverM(x[0]);
234 Double_t mass = par[0];
235 Int_t charge = par[1];
236 if (mass < 0) {mass = - mass; scale = 2;}
237 Double_t poverm = pove/mass;
241 dx2 = TMath::Log2(5.);
243 scale *= charge*charge;
244 return TMath::Log10(scale*TMath::Exp(Bichsel::Instance()->I70Trs(part,TMath::Log10(poverm))));
247 Double_t bichselZTrs(Double_t *x,Double_t *par) {
248 Double_t pove = pOverM(x[0]);
250 Double_t mass = par[0];
252 if (part < 0 || part > 9) part = 0;
253 if (mass < 0) {mass = - mass; scale = 2;}
254 Double_t poverm = pove/mass;
255 Double_t charge = 1.;
260 dx2 = TMath::Log2(5.);
262 scale *= charge*charge;
263 return TMath::Log10(scale*TMath::Exp(Bichsel::Instance()->IfitTrs(part,TMath::Log10(poverm))));
267 Double_t dNdx(Double_t *x,Double_t *par) {
268 Double_t pove = pOverM(x[0]);
270 Double_t mass = par[0];
271 if (mass < 0) {mass = - mass; scale = 2;}
272 Double_t poverm = pove/mass;
273 Double_t charge = 1.;
275 if (par[1] > 1.0) charge = par[1];
277 return TMath::Log10(scale*StdEdxPull::EvalPred(poverm, 2, charge));
280 Double_t dNdxOld(Double_t *x,Double_t *par) {
281 Double_t pove = pOverM(x[0]);
283 Double_t mass = par[0];
284 if (mass < 0) {mass = - mass; scale = 2;}
285 Double_t poverm = pove/mass;
286 Double_t charge = 1.;
288 if (par[1] > 1.0) charge = par[1];
290 return TMath::Log10(scale*StdEdxPull::EvalPred2(poverm, dx2, 2, charge));
292 #if !defined(__CINT__) && !defined(__CLING__)
294 Double_t aleph70P(Double_t *x,Double_t *par) {
301 Double_t b2inv = 1. + 1./(bg*bg);
302 Double_t beta = 1./TMath::Sqrt(b2inv);
303 Double_t dEdx = par[0]*(-1 + TMath::Power(beta,-par[3])*(par[1] - TMath::Log(TMath::Max(1e-10,par[2] + TMath::Power(bg,-par[4])))));
307 Double_t aleph70(Double_t *x,Double_t *par) {
308 static const Double_t dEdxMIP = 2.39761562607903311;
309 static Double_t MIPBetaGamma = 4.;
311 struct Par_t {Int_t h, N; Double_t xmin, xmax, pars[10];};
312 const Par_t Par[9] = {
314 { 0, 4, 3.0, 6.0,{ 0.14105, -0.09078, 0.01901, -0.00128, 0, 0, 0, 0, 0, 0}},
315 { 1, 2, 0.0, 4.5,{ -0.00689, 0.00256, 0, 0, 0, 0, 0, 0, 0, 0}},
316 { 2, -1, 0.0, 4.5,{ 0.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
317 { 3, 9, -0.1, 3.7,{ 0.00869, 0.21918, -0.88919, 1.30023, -0.97075, 0.41298, -0.10214, 0.01379, -0.00079, 0}},
318 { 4, 9, -0.6, 3.3,{ 0.03052, -0.02423, -0.05636, -0.11585, 0.41292, -0.38956, 0.16837, -0.03494, 0.00283, 0}},
319 { 5, 8, -1.0, 2.9,{ 0.03523, -0.10625, 0.04182, 0.07820, -0.03816, -0.02735, 0.01940, -0.00304, 0, 0}},
320 { 6, 10, -1.0, 2.8,{ 0.03092, -0.07846, 0.01668, 0.00331, 0.12771, -0.05480, -0.13897, 0.13928, -0.04715, 0.00555}},
321 { 7, 9, -0.8, 2.9,{ 0.09158, -0.07816, 0.07039, 0.00578, -0.16160, 0.26547, -0.18728, 0.05988, -0.00710, 0}},
322 { 8, 10, -0.8, 2.9,{ 0.09366, -0.08276, 0.06191, 0.02631, -0.17044, 0.30536, -0.26867, 0.11847, -0.02505, 0.00201}}
325 static Double_t ppar[7] = { 0.0857988, 9.46138, 0.000206655, 2.12222, 0.974, -1, 0.13957};
329 static Double_t ppar[7] = {0.12337, 6.61371, 0.00201, 2.27381, 1.54174, -1.00000, 0 };
331 static Double_t Norm = dEdxMIP/aleph70P(&MIPBetaGamma,ppar);
332 Int_t hyp = (Int_t ) par[0];
333 Int_t h = Part[hyp].Index;
334 Double_t ScaleL10 = 0;
337 ScaleL10 = TMath::Log10(2.);
339 Double_t pove = pOverM(x[0]);
340 Double_t mass = Part[hyp].Mass;
341 Double_t poverm = pove/mass;
342 Double_t charge = 1.;
343 if (h > 6 && h > 9) charge = 2;
344 else if (h == 10 || h == 11) charge = 3;
346 Double_t bg = poverm;
352 Double_t dEdxL10 = TMath::Log10(Norm*aleph70P(&bg,ppar));
355 TString fName(Form(
"pol%i",Par[h].N-1));
356 TF1 *f = (TF1 *) gROOT->GetListOfFunctions()->FindObject(fName);
358 f =
new TF1(fName,fName,0,1);
360 f->SetParameters(&Par[h].pars[0]);
361 Double_t bgL10 = TMath::Log10(bg);
362 dEdxL10 += f->Eval(bgL10);
365 return 2*TMath::Log10(charge) + dEdxL10 + ScaleL10;
369 void bichselG10(
const Char_t *type=
"zN", Int_t Nhyps = 9, Bool_t rigidity = kFALSE) {
370 fgRigidity = rigidity;
371 if (gClassTable->GetID(
"StBichsel") < 0) {
372 gSystem->Load(
"libTable");
373 gSystem->Load(
"St_base");
374 gSystem->Load(
"StarClassLibrary");
375 gSystem->Load(
"StBichsel");
379 TLegend *leg =
new TLegend(0.85,0.45,0.95,0.9,
"");
382 for (Int_t i = NF-1; i >=0; i--) {
383 if (Type.Contains(FNames[i],TString::kIgnoreCase)) {
388 if (Nhyps < 9) Nhyps = 9;
389 if (Nhyps > NMasses) Nhyps = NMasses;
390 TH1D *FHyps[20] = {0};
391 for (
int h = 0; h < Nhyps; h++) {
394 Char_t *FunName = Form(
"%s%s%i",FNames[f],Part[h].Name,(
int)log2dx[dx]);
395 cout <<
"Make " << h <<
"\t" << FunName << endl;
396 Double_t xmin = -xmax;
401 if (f == 3) func =
new TF1(FunName,bichsel70,xmin, xmax,3);
402 else if (f == 2) func =
new TF1(FunName,bichselZ ,xmin, xmax,3);
403 else if (f == 5) func =
new TF1(FunName,bichsel70M ,xmin, xmax,3);
404 else if (f == 6) func =
new TF1(FunName,dNdx ,xmin, xmax,3);
405 else if (f == 7) func =
new TF1(FunName,bichselZM,xmin, xmax,3);
406 else if (f == 8) func =
new TF1(FunName,dEdxModelZ,xmin, xmax,3);
408 else if (f == 9) func =
new TF1(FunName,bichsel70Trs,xmin, xmax,3);
409 else if (f ==10) func =
new TF1(FunName,bichselZTrs,xmin, xmax,3);
414 func->SetParameter(0,Part[h].Mass);
415 func->SetParameter(1,Part[h].Charge);
417 if (color > 8) color -= 8;
420 func->SetLineColor(color);
421 func->SetMarkerColor(color);
424 leg->AddEntry(func,Part[h].Name);
425 #if !defined( __CINT__) && defined(__Aleph__)
426 TF1 *fA =
new TF1(Form(
"Aleph%s",Part[h].Name),aleph70,xmin,xmax, 1);
427 fA->SetParameter(0,h);
428 fA->SetLineColor(color);
429 fA->SetMarkerColor(color);
432 leg->AddEntry(fA,Part[h].Name);
434 FHyps[h] = (TH1D *) func->GetHistogram();
435 FHyps[h]->SetName(func->GetName());
439 TH1D *FHypsFromPion[20] = {0};
440 if (FHyps[kpionIdx]) {
441 TCanvas *csep = (TCanvas *) gROOT->GetListOfCanvases()->FindObject(Form(
"c2",type));
442 if (! csep) csep =
new TCanvas(Form(
"c2",type),Form(
"Separation %s",type));
444 TLegend *legS =
new TLegend(0.85,0.45,0.95,0.9,
"");
445 TH1F *frame = csep->DrawFrame(-1,-1,4,2);
447 for (
int h = 0; h < Nhyps; h++) {
448 if (! FHyps[h] )
continue;
449 FHypsFromPion[h] =
new TH1D(*FHyps[h]);
450 FHypsFromPion[h]->SetName(Form(
"%s-%s",FHyps[h]->GetName(),FHyps[kpionIdx]->GetName()));
451 FHypsFromPion[h]->Add(FHyps[kpionIdx],-1);
452 if (Type.Contains(
"dNdx",TString::kIgnoreCase)) FHypsFromPion[h]->Scale(1.1);
453 FHypsFromPion[h]->Draw(
"samel");
454 legS->AddEntry(FHypsFromPion[h],FHypsFromPion[h]->GetName());