StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
bichselG10.C
1 /*
2  root.exe lBichsel.C bichselG10.C+
3  bichselG10("N"); // dN/dx
4  bichselG10("70"; // I70
5  bichselG10("z"); // Ifit
6  TH1D *pB70 = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("70p1"))->GetHistogram();
7  TH1D *piB70 = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("70#pi1"))->GetHistogram();
8  TH1D *diffB70 = new TH1D(*pB70);
9  diffB70->SetName("diffB70");
10  diffB70->Add(pB70,piB70,1,-1);
11  diffB70->SetLineColor(1);
12  TH1D *pBz = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("zp1"))->GetHistogram();
13  TH1D *piBz = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("z#pi1"))->GetHistogram();
14  TH1D *diffBz = new TH1D(*pBz);
15  diffBz->SetName("diffBz");
16  diffBz->Add(pBz,piBz,1,-1);
17  diffBz->SetLineColor(2);
18  TH1D *pdNdx = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("dNdxp1"))->GetHistogram();
19  TH1D *pidNdx = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("dNdx#pi1"))->GetHistogram();
20  TH1D *diffdNdx = new TH1D(*pdNdx);
21  diffdNdx->SetName("diffdNdx");
22  diffdNdx->Add(pdNdx,pidNdx,1,-1);
23  diffdNdx->SetLineColor(3);
24  cppi = new TCanvas("cppi","cppi");
25  diffB70->SetXTitle("log_{10}p")
26  diffB70->SetTitle("z_{p} - z_{#pi}");
27  diffB70->Draw("l");
28  diffBz->Draw("samel");
29  diffdNdx->Draw("samel");
30  TLegend *l = new TLegend(0.6,0.6,0.8,0.8);
31  l->AddEntry(diff70,"I70");
32  l->AddEntry(diffBz,"Ifit");
33  l->AddEntry(diffdNdx,"dNdx");
34 
35  TH1D *eB70 = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("70e1"))->GetHistogram();
36  TH1D *piB70 = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("70#pi1"))->GetHistogram();
37  TH1D *diffB70 = new TH1D(*eB70);
38  diffB70->SetName("diffB70");
39  diffB70->Add(eB70,piB70,1,-1);
40  diffB70->SetLineColor(1);
41  TH1D *eBz = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("ze1"))->GetHistogram();
42  TH1D *piBz = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("z#pi1"))->GetHistogram();
43  TH1D *diffBz = new TH1D(*eBz);
44  diffBz->SetName("diffBz");
45  diffBz->Add(eBz,piBz,1,-1);
46  diffBz->SetLineColor(2);
47  TH1D *edNdx = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("dNdxe1"))->GetHistogram();
48  TH1D *pidNdx = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("dNdx#pi1"))->GetHistogram();
49  TH1D *diffdNdx = new TH1D(*edNdx);
50  diffdNdx->SetName("diffdNdx");
51  diffdNdx->Add(edNdx,pidNdx,1,-1);
52  diffdNdx->SetLineColor(3);
53  cepi = new TCanvas("cepi","cepi");
54  diffB70->SetXTitle("log_{10}p")
55  diffB70->SetTitle("z_{e} - z_{#pi}");
56  diffB70->Draw("l");
57  diffBz->Draw("samel");
58  diffdNdx->Draw("samel");
59  TLegend *l = new TLegend(0.6,0.6,0.8,0.8);
60  l->AddEntry(diff70,"I70");
61  l->AddEntry(diffBz,"Ifit");
62  l->AddEntry(diffdNdx,"dNdx");
63 
64 */
65 #if !defined(__CINT__)
66 // code that should be seen ONLY by the compiler
67 #else
68 #if !defined(__CINT__) || defined(__MAKECINT__)
69 // code that should be seen by the compiler AND rootcint
70 #else
71 // code that should always be seen
72 #endif
73 #endif
74 //________________________________________________________________________________
75 #if !defined(__CINT__) || defined(__MAKECINT__)
76 #include "Riostream.h"
77 #include <stdio.h>
78 #include "TF1.h"
79 #include "TMath.h"
80 #include "TSystem.h"
81 #include "TCanvas.h"
82 #include "TClassTable.h"
83 #include "StBichsel/Bichsel.h"
84 #include "StBichsel/StdEdxModel.h"
85 #include "StBichsel/StdEdxPull.h"
86 #include "TLegend.h"
87 #include "TROOT.h"
88 #include "TString.h"
89 #else
90 class Bichsel;
91 #endif
92 // function object (functor)
93 struct MyDerivFunc {
94  MyDerivFunc(TF1 * f): fFunc(f) {}
95  double operator() (double *x, double * ) const {
96  return fFunc->Derivative(*x);
97  }
98  TF1 * fFunc;
99 };
100 //________________________________________________________________________________
101 TF1 *Derivative(TF1 *f1) {
102  MyDerivFunc * deriv = new MyDerivFunc(f1);
103  TString name("der_");
104  name += f1->GetName();
105  TF1 * f2 = new TF1(name,deriv, f1->GetXmin(), f1->GetXmax(), 0, "MyDerivFunc");
106  return f2;
107 }
108 //Bichsel *m_Bichsel = 0;
109 const Int_t NMasses = 20;
110 const Double_t kAu2Gev=0.9314943228;
111 struct Part_t {
112  const Char_t *Name;
113  Int_t PiD;
114  Int_t Charge;
115  Int_t Index;
116  Double_t Mass;
117 };
118 const Int_t kpionIdx = 2;
119 Part_t Part[NMasses] = {// https://periodictable.com/Isotopes/
120  //name, PiD, Charge, Index, Mass
121  {"e", 0, 1, 0, 0.51099907e-3}, // 3 e
122  {"#mu", 4, 1, 1, 0.1056584}, // 5 mu
123  {"#pi", 3, 1, 2, 0.13956995}, // 2 pi
124  {"K", 2, 1, 3, 0.493677}, // 1 K
125  {"p", 1, 1, 4, 0.93827231}, // 0 p
126  {"d", 5, 1, 5, 1.87561339}, // 4 d
127  {"t", 6, 1, 6, 2.80925}, // 6 t
128  {"He3", 7, 2, 7, 2.80923}, //GEANT3 // 7 He3
129  {"#alpha", 8, 2, 8, 4.00260325415*kAu2Gev}, // 8 He4
130  {"He6", -1, 2, 12, 6.018889124*kAu2Gev}, // 9
131  {"He8", -1, 2, 12, 8.033921897*kAu2Gev}, //10
132  {"Li6", 10, 3, 10, 5.6031}, //11 Li6
133  {"Li7", 9, 3, 9, 6.5354}, //12 Li7
134  {"Be7", -1, 4, 0, 7.016003437*kAu2Gev}, //13
135  {"Be9", -1, 4, 0, 9.01218307*kAu2Gev}, //14
136  {"Be10", -1, 4, 0, 10.01353470*kAu2Gev}, //15
137  {"B10", -1, 5, 0, 10.01353470*kAu2Gev}, //16
138  {"B11", -1, 5, 0, 11.009305167*kAu2Gev}, //17
139  {"2#pi", 3, -1, -2, -0.13956995}, //18 2*pi
140  {"2p", 1, -1, 0, -0.93827231} //19 2*p
141 };
142 const Int_t NF = 9; // 0, 1, 2, 3, 4, 5. 6, 7, 8. 9, 10,
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;
147 //________________________________________________________________________________
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);
151 }
152 //________________________________________________________________________________
153 Double_t bichselZ(Double_t *x,Double_t *par) {
154  Double_t pove = pOverM(x[0]);
155  Double_t scale = 1;
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.;
160  Double_t dx2 = 1;
161  if (par[1] > 1.0) {
162  charge = par[1];
163  poverm *= charge;
164  dx2 = TMath::Log2(5.);
165  }
166  return TMath::Log10(scale*charge*charge*TMath::Exp(Bichsel::Instance()->GetMostProbableZ(TMath::Log10(poverm),dx2)));//TMath::Exp(7.81779499999999961e-01));
167  // Charge*Charge* (TMath::Exp(Bichsel::Instance()->GetMostProbableZM(TMath::Log10(TMath::Abs(Charge)*p/M),dx2)))
168  // return TMath::Log10(1e6*scale*StdEdxPull::EvalPred(poverm, 1, charge));
169 }
170 //________________________________________________________________________________
171 Double_t dEdxModelZ(Double_t *x,Double_t *par) { //new dEdxModel
172  Double_t pove = pOverM(x[0]);
173  Double_t scale = 1;
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));
180 }
181 //________________________________________________________________________________
182 Double_t bichselZM(Double_t *x,Double_t *par) {
183  Double_t pove = pOverM(x[0]);
184  Double_t scale = 1;
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.;
189  Double_t dx2 = 1;
190  if (par[1] > 1.0) {
191  charge = par[1];
192  poverm *= TMath::Abs(charge);
193  dx2 = TMath::Log2(5.);
194  }
195  return TMath::Log10(1e6*scale*StdEdxPull::EvalPred(poverm, 1, charge));
196 }
197 //________________________________________________________________________________
198 Double_t bichsel70(Double_t *x,Double_t *par) {
199  Double_t pove = pOverM(x[0]);
200  Double_t scale = 1;
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.;
205  Double_t dx2 = 1;
206  if (par[1] > 1.0) {
207  charge = par[1];
208  poverm *= charge;
209  dx2 = TMath::Log2(5.);
210  }
211  return TMath::Log10(1e6*scale*StdEdxPull::EvalPred(poverm, 0, charge));
212 }
213 //________________________________________________________________________________
214 Double_t bichsel70M(Double_t *x,Double_t *par) {
215  Double_t pove = pOverM(x[0]);
216  Double_t scale = 1;
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.;
221  Double_t dx2 = 1;
222  if (par[1] > 1.0) {
223  charge = par[1];
224  poverm *= charge;
225  dx2 = TMath::Log2(5.);
226  }
227  return TMath::Log10(1e6*scale*StdEdxPull::EvalPred(poverm, 0, charge));
228 }
229 #if 0
230 //________________________________________________________________________________
231 Double_t bichsel70Trs(Double_t *x,Double_t *par) {
232  Double_t pove = pOverM(x[0]);
233  Double_t scale = 1;
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;
238  Double_t dx2 = 1;
239  if (charge > 1) {
240  poverm *= charge;
241  dx2 = TMath::Log2(5.);
242  }
243  scale *= charge*charge;
244  return TMath::Log10(scale*TMath::Exp(Bichsel::Instance()->I70Trs(part,TMath::Log10(poverm))));
245 }
246 //________________________________________________________________________________
247 Double_t bichselZTrs(Double_t *x,Double_t *par) {
248  Double_t pove = pOverM(x[0]);
249  Double_t scale = 1;
250  Double_t mass = par[0];
251  Int_t part = par[2];
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.;
256  Double_t dx2 = 1;
257  if (par[1] > 1.0) {
258  charge = par[1];
259  poverm *= charge;
260  dx2 = TMath::Log2(5.);
261  }
262  scale *= charge*charge;
263  return TMath::Log10(scale*TMath::Exp(Bichsel::Instance()->IfitTrs(part,TMath::Log10(poverm))));//TMath::Exp(7.81779499999999961e-01));
264 }
265 #endif
266 //________________________________________________________________________________
267 Double_t dNdx(Double_t *x,Double_t *par) {
268  Double_t pove = pOverM(x[0]);
269  Double_t scale = 1;
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.;
274  Double_t dx2 = 1;
275  if (par[1] > 1.0) charge = par[1];
276  poverm *= charge;
277  return TMath::Log10(scale*StdEdxPull::EvalPred(poverm, 2, charge));
278 }
279 //________________________________________________________________________________
280 Double_t dNdxOld(Double_t *x,Double_t *par) {
281  Double_t pove = pOverM(x[0]);
282  Double_t scale = 1;
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.;
287  Double_t dx2 = 1;
288  if (par[1] > 1.0) charge = par[1];
289  poverm *= charge;
290  return TMath::Log10(scale*StdEdxPull::EvalPred2(poverm, dx2, 2, charge));
291 }
292 #if !defined(__CINT__) && !defined(__CLING__)
293 //________________________________________________________________________________
294 Double_t aleph70P(Double_t *x,Double_t *par) {
295  /*
296  W.Blum, L. Rolandi "Particle Detection with Drift Chambers", page 246, eq. (9.5)
297  F_g(v) = p[0]/beta**p[3]*(p[1] - beta**p[3] - log(p[2] + (beta*gamma)**-p[4]);
298  F_g(v) = p[0]*(1/beta**p[3]*(p[1] - log(p[2] + 1/(beta*gamma)**p[4])) - 1)
299  */
300  Double_t bg = x[0];
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])))));
304  return dEdx;
305 };
306 //________________________________________________________________________________
307 Double_t aleph70(Double_t *x,Double_t *par) {
308  static const Double_t dEdxMIP = 2.39761562607903311;
309  static Double_t MIPBetaGamma = 4.;
310 #if 0
311  struct Par_t {Int_t h, N; Double_t xmin, xmax, pars[10];};
312  const Par_t Par[9] = {
313  /* name h n+1 xmin xmax pars[n+1] */
314  /* electron */{ 0, 4, 3.0, 6.0,{ 0.14105, -0.09078, 0.01901, -0.00128, 0, 0, 0, 0, 0, 0}},
315  /* muon */{ 1, 2, 0.0, 4.5,{ -0.00689, 0.00256, 0, 0, 0, 0, 0, 0, 0, 0}},
316  /* pion */{ 2, -1, 0.0, 4.5,{ 0.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
317  /* kaon */{ 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  /* proton */{ 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  /* deuteron */{ 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  /* triton */{ 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  /* He3 */{ 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  /* alpha */{ 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}}
323  };
324  // const Double_t ppar[7] = { 0.0857988, 9.46138, 0.000206655, 2.12222, 0.974, -1, 0.13957}; /* pion */
325  static Double_t ppar[7] = { 0.0857988, 9.46138, 0.000206655, 2.12222, 0.974, -1, 0.13957}; /* pion */
326 #else
327  // static Double_t ppar[7] = { 0.0762, 10.632, 0.134e-4, 1.863, 1.948, -1, -1}; /* Aleph parameters from Alice TPC TDR */
328  // static Double_t ppar[7] = { 0.08942, 8.91971, 0.00024, 2.27383, 1.54174, -1.00000, 0}; /* pion */
329  static Double_t ppar[7] = {0.12337, 6.61371, 0.00201, 2.27381, 1.54174, -1.00000, 0 }; /* g All */
330 #endif
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;
335  if (h < 0) {
336  h = -h;
337  ScaleL10 = TMath::Log10(2.);
338  }
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;
345  poverm *= charge;
346  Double_t bg = poverm;
347  /*
348  W.Blum, L. Rolandi "Particle Detection with Drift Chambers", page 246, eq. (9.5)
349  F_g(v) = p[0]/beta**p[3]*(p[1] - beta**p[3] - log(p[2] + (beta*gamma)**-p[4]));
350  F_g(v) = p[0]*(1/beta**p[3]*(p[1] - log(p[2] + 1/(beta*gamma)**p[4])) - 1)
351  */
352  Double_t dEdxL10 = TMath::Log10(Norm*aleph70P(&bg,ppar));
353 #if 0
354  if (Par[h].N > 0) {
355  TString fName(Form("pol%i",Par[h].N-1));
356  TF1 *f = (TF1 *) gROOT->GetListOfFunctions()->FindObject(fName);
357  if (! f) {
358  f = new TF1(fName,fName,0,1);
359  }
360  f->SetParameters(&Par[h].pars[0]);
361  Double_t bgL10 = TMath::Log10(bg);
362  dEdxL10 += f->Eval(bgL10);
363  }
364 #endif
365  return 2*TMath::Log10(charge) + dEdxL10 + ScaleL10;
366 }
367 #endif /* __CINT__ */
368 //________________________________________________________________________________
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");
376  }
377  // m_Bichsel = Bichsel::Instance();
378  TString Type(type);
379  TLegend *leg = new TLegend(0.85,0.45,0.95,0.9,"");
380  Double_t xmax = 4;
381  Int_t f = 3;
382  for (Int_t i = NF-1; i >=0; i--) {
383  if (Type.Contains(FNames[i],TString::kIgnoreCase)) {
384  f = i;
385  break;
386  }
387  }
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++) { // Masses
392  // for (int h = 0; h < 7; h++) { // Masses
393  Int_t dx = 1;
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;
397  // if (h == 0 || h >= 5) xmin = -0.75;
398  // if (h == 4) xmin = -0.70;
399  // if (h == 6) xmin = -0.50;
400  TF1 *func = 0;
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);
407 #if 0
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);
410 #endif
411  else {
412  return;
413  }
414  func->SetParameter(0,Part[h].Mass);
415  func->SetParameter(1,Part[h].Charge);
416  Int_t color = h+1;
417  if (color > 8) color -= 8;
418  // if (color > 7) color++;
419 #if 1
420  func->SetLineColor(color);
421  func->SetMarkerColor(color);
422 #endif
423  func->Draw("same");
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);
430  fA->SetLineStyle(2);
431  fA->Draw("same");
432  leg->AddEntry(fA,Part[h].Name);
433 #endif
434  FHyps[h] = (TH1D *) func->GetHistogram();
435  FHyps[h]->SetName(func->GetName());
436  }
437  leg->Draw();
438 #if 0
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));
443  else csep->Clear();
444  TLegend *legS = new TLegend(0.85,0.45,0.95,0.9,"");
445  TH1F *frame = csep->DrawFrame(-1,-1,4,2);
446  // separation from pion
447  for (int h = 0; h < Nhyps; h++) { // Masses
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());
455  }
456  legS->Draw();
457  csep->Update();
458  }
459 #endif
460 }