All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
LikelihoodPID.cc
Go to the documentation of this file.
1 /*
2 --------------- inplement Particle Identofication using Global likelihood -------------
3 
4 This code includes class for Particle ID. The method is Naive Bayes Classfier
5 1. First, check whether a particle is electron or not.
6  If so, the particle is identified as electron and return electron status code.
7 2. And then, check whether a particle a particle is muonn or not.
8  If so, the particle is identified as muon and return muon status code.
9 3. If a particle is neither electron nor muon,
10  try to classify a particle to each Hadron type(Pion, Kaon, or Proton)
11 
12 TODO:
13 Bethe-Bloch parameters should be moved to steer file
14 Threshold parameters should be moved to steer file
15 
16 make flags to choose Particle ID method
17 
18 risk-minimization and MAP
19  */
20 
21 #include <string>
22 #include <sstream>
23 #include <TFile.h>
24 #include <TH1F.h>
25 
26 #include "DD4hep/Detector.h"
27 #include "DDRec/Vector3D.h"
28 #include "DD4hep/DD4hepUnits.h"
29 
30 #include "LikelihoodPID.hh"
31 
32 using namespace std;
33 
34 std::string itos(int i)
35 {
36  std::stringstream s;
37  s << i;
38  return s.str();
39 }
40 
41 LikelihoodPID::LikelihoodPID(double *pars){
42  //set parameters for bethebloch
43  //electron
44  for(int i=0;i<5;i++) par[0][i]=pars[i];
45  //muon
46  for(int i=0;i<5;i++) par[1][i]=pars[i+5];
47  //pion
48  for(int i=0;i<5;i++) par[2][i]=pars[i+10];
49  //kaon
50  for(int i=0;i<5;i++) par[3][i]=pars[i+15];
51  //proton
52  for(int i=0;i<5;i++) par[4][i]=pars[i+20];
53 
54  //choose method
55  _usebayes=(int)pars[25];
56  _dEdxnorm=(float)pars[26];
57  _dEdxerrfact=pars[27];
58 
59  //set mass
60  emass=0.000510998;
61  mmass=0.105658;
62  pimass=0.139570;
63  kmass=0.493677;
64  pmass=0.938272;
65 
66  //get B field
67  double bfield[3]={};
68  dd4hep::Detector& lcdd = dd4hep::Detector::getInstance();
69  lcdd.field().magneticField({0.0,0.0,0.0}, bfield);
70  _bfield = (float)bfield[2]/dd4hep::tesla;
71 
72  return;
73 }
74 
75 LikelihoodPID::LikelihoodPID(string fname, double *pars, std::vector<float> cost){
76  //set parameters for bethebloch
77  //electron
78  for(int i=0;i<5;i++) par[0][i]=pars[i];
79  //muon
80  for(int i=0;i<5;i++) par[1][i]=pars[i+5];
81  //pion
82  for(int i=0;i<5;i++) par[2][i]=pars[i+10];
83  //kaon
84  for(int i=0;i<5;i++) par[3][i]=pars[i+15];
85  //proton
86  for(int i=0;i<5;i++) par[4][i]=pars[i+20];
87 
88  //choose method
89  _usebayes=(int)pars[25];
90  _dEdxnorm=(float)pars[26];
91  _dEdxerrfact=pars[27];
92 
93  //set mass
94  emass=0.000510998;
95  mmass=0.105658;
96  pimass=0.139570;
97  kmass=0.493677;
98  pmass=0.938272;
99 
100  //get pdf plots
101  string ffstr1 = fname; // "./pdf_standard_s_05_v01.root";
102  fpdf=new TFile(ffstr1.c_str());
103 
104  string hname,hname2;
105  //for(int i=0;i<6;i++) pdf[i] =new TH1F()[14];
106  for(int i=0;i<6;i++){
107  //ep
108  hname="hep" + itos(i+1);
109  hname2="hep" + itos(i+1) + "_2";
110  pdf[i][0]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
111  //ehad
112  hname="hehad" + itos(i+1);
113  hname2="hehad" + itos(i+1) + "_2";
114  pdf[i][1]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
115  //mucal
116  hname="hmucal" + itos(i+1);
117  hname2="hmucal" + itos(i+1) + "_2";
118  pdf[i][2]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
119  //chi2
120  hname="hchi2" + itos(i+1);
121  hname2="hchi2" + itos(i+1) + "_2";
122  pdf[i][3]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
123  //showermax/exp.showermax
124  hname="hldiscrepancy" + itos(i+1);
125  hname2="hldiscrepancy" + itos(i+1) + "_2";
126  pdf[i][4]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
127  //absorption length
128  hname="htdiscrepancy" + itos(i+1);
129  hname2="htdiscrepancy" + itos(i+1) + "_2";
130  pdf[i][5]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
131  //xl20
132  hname="hxl20" + itos(i+1);
133  hname2="hxl20" + itos(i+1) + "_2";
134  pdf[i][6]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
135  //likeliele
136  hname="hlikeliele" + itos(i+1);
137  hname2="hlikeliele" + itos(i+1) + "_2";
138  pdf[i][7]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
139  //likelimuo
140  hname="hlikelimuo" + itos(i+1);
141  hname2="hlikelimuo" + itos(i+1) + "_2";
142  pdf[i][8]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
143  //likelipi
144  hname="hlikelipi" + itos(i+1);
145  hname2="hlikelipi" + itos(i+1) + "_2";
146  pdf[i][9]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
147  //likelik
148  hname="hlikelik" + itos(i+1);
149  hname2="hlikelik" + itos(i+1) + "_2";
150  pdf[i][10]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
151  //likelip
152  hname="hlikelip" + itos(i+1);
153  hname2="hlikelip" + itos(i+1) + "_2";
154  pdf[i][11]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
155 
156  TObject *obj;
157  //cluster depth
158  hname="hcludepth" + itos(i+1);
159  fpdf->GetObject(hname.c_str(), obj);
160  if(obj){
161  hname2="hcludepth" + itos(i+1) + "_2";
162  pdf[i][14]=(TH1F*)obj->Clone(hname2.c_str());
163 
164  //deltax
165  hname="hdeltax" + itos(i+1);
166  fpdf->GetObject(hname.c_str(), obj);
167  hname2="hdeltax" + itos(i+1) + "_2";
168  pdf[i][12]=(TH1F*)obj->Clone(hname2.c_str());
169 
170  //deltaz
171  hname="hdeltaz" + itos(i+1);
172  fpdf->GetObject(hname.c_str(), obj);
173  hname2="hdeltaz" + itos(i+1) + "_2";
174  pdf[i][13]=(TH1F*)obj->Clone(hname2.c_str());
175 
176  //hitmean
177  hname="hhitmean" + itos(i+1);
178  fpdf->GetObject(hname.c_str(), obj);
179  hname2="hhitmean" + itos(i+1) + "_2";
180  pdf[i][15]=(TH1F*)obj->Clone(hname2.c_str());
181 
182  //hitrms
183  hname="hhitrms" + itos(i+1);
184  fpdf->GetObject(hname.c_str(), obj);
185  hname2="hhitrms" + itos(i+1) + "_2";
186  pdf[i][16]=(TH1F*)obj->Clone(hname2.c_str());
187  }else{
188  pdf[i][12]=NULL; //backward compatibility
189  pdf[i][13]=NULL; //backward compatibility
190  pdf[i][14]=NULL; //backward compatibility
191  pdf[i][15]=NULL; //backward compatibility
192  pdf[i][16]=NULL; //backward compatibility
193  }
194 
195  //deltaz
196  /*hname="hdeltaz" + itos(i+1);
197  hname2="hdeltaz" + itos(i+1) + "_2";
198  pdf[i][13]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
199  //cluster depth
200  hname="hcludepth" + itos(i+1);
201  hname2="hcludepth" + itos(i+1) + "_2";
202  pdf[i][14]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
203  //hitmean
204  hname="hhitmean" + itos(i+1);
205  hname2="hhitmean" + itos(i+1) + "_2";
206  pdf[i][15]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());
207  //hitrms
208  hname="hhitrms" + itos(i+1);
209  hname2="hhitrms" + itos(i+1) + "_2";
210  pdf[i][16]=(TH1F*)fpdf->Get(hname.c_str())->Clone(hname2.c_str());*/
211  }
212 
213  //normalize histograms
214  double weight=1.0;
215  for(int i=0;i<6;i++){
216  for(int j=0;j<17;j++){
217  //normalize histograms
218  if(pdf[i][j]!=NULL){
219  weight=pdf[i][j]->Integral(0,pdf[i][j]->GetNbinsX()+1,"");
220  pdf[i][j]->Scale(1.0/weight);
221  _weights[i][j]=weight; //for hadron likelihood calculation
222  }
223  }
224  }
225 
226  //normalize weights
227  for(int j=0;j<17;j++){
228  if(pdf[2][j]!=NULL && pdf[3][j]!=NULL && pdf[4][j]!=NULL){
229  double denom = _weights[2][j]+_weights[3][j]+_weights[4][j];
230  for(int i=2;i<5;i++) _weights[i][j] = _weights[i][j]/denom;
231  }else
232  for(int i=2;i<5;i++) _weights[i][j] = 1.0 / 3.0; //same weights
233  }
234 
235  //set threshold
236  //this is original
237  threshold[0]=0.0; //TMath::Exp(-0.55);
238  threshold[1]=0.0; //TMath::Exp(-0.7);
239  threshold[2]=0.0;
240  threshold[3]=0.0; //TMath::Exp(-0.6);
241  threshold[4]=0.0; //TMath::Exp(-0.6);
242 
243  //set cost
244  if(_usebayes==1){
245  for(int i=0;i<5;i++){
246  for(int j=0;j<5;j++){
247  if(i==j) penalty[i][j]=1.0e-50;
248  else penalty[i][j]=1.0;
249  }
250  }
251  }else if(_usebayes==2){
252  for(int i=0;i<5;i++){
253  for(int j=0;j<5;j++){
254  penalty[i][j]=cost[i*5+j];
255  }
256  }
257  }
258 
259  //get z component of B field
260  double bfield[3]={};
261  dd4hep::Detector& lcdd = dd4hep::Detector::getInstance();
262  lcdd.field().magneticField({0.0,0.0,0.0}, bfield);
263  _bfield = (float)bfield[2]/dd4hep::tesla;
264  //cout << "magnetic field: " << _bfield << endl;
265 
266  return;
267 }
268 
269 LikelihoodPID::~LikelihoodPID(){
270 
271  fpdf->Close();
272 
273  return;
274 }
275 
276 //public
277 int LikelihoodPID::Classification(TLorentzVector pp, EVENT::Track* trk, EVENT::ClusterVec& cluvec){
278  int tmpid=-1;
279 
280  //set prior first
281  for(int i=0;i<5;i++){
282  prior[i]=0.2;
283  _posterior[i]=prior[i];
284  _likelihood[i]=0.0;
285  }
286  _posterior[5]=prior[2]+prior[3]+prior[4];
287  _likelihood[5]=0.0;
288 
289  //get dedx distance
290  for(int i=0;i<5;i++) _dEdxDist[i]=get_dEdxChi2(i, pp.Vect(), trk->getdEdxError(),trk->getdEdx());
291 
292  //first check whether electron or not
293  for(int i=0;i<5;i++){
294  _likelihood[i]=0.0;
295  if(_basicFlg==false && _dEdxFlg==true && _showerShapesFlg==false) prior[i]=0.2;
296  }
297  _likelihood[5]=0.0; //hadron type
298 
299  tmpid=Class_electron(pp, trk, cluvec);
300  if(tmpid==0) return tmpid;
301 
302  //second check whether muon or not
303  for(int i=0;i<5;i++){
304  _likelihood[i]=0.0;
305  if(_basicFlg==false && _dEdxFlg==true && _showerShapesFlg==false) prior[i]=0.2;
306  }
307  _likelihood[5]=0.0; //hadron type
308 
309  tmpid=Class_muon(pp, trk, cluvec);
310  if(tmpid==1) return tmpid;
311 
312  //third classify hadrons
313  for(int i=0;i<5;i++){
314  _likelihood[i]=0.0;
315  if(_basicFlg==false && _dEdxFlg==true && _showerShapesFlg==false) prior[i]=0.2;
316  }
317  _likelihood[5]=0.0; //hadron type
318 
319  tmpid=Class_hadron(pp, trk, cluvec);
320 
321  //avoid strange value
322  if(_posterior[0]!=_posterior[0]){ //cannot estimate likelihood and probability
323  for(int i=0;i<6;i++){
324  _likelihood[i] = 999.0;
325  _posterior[i] = 0.0;
326  }
327  }else if(_likelihood[0]==0.0 &&_likelihood[1]==0.0 &&_likelihood[2]==0.0 && _likelihood[3]==0.0 && _likelihood[4]==0.0){
328  for(int i=0;i<6;i++){
329  _likelihood[i] = 999.0;
330  _posterior[i] = 0.0;
331  }
332  }
333 
334  return tmpid;
335 }
336 
337 double *LikelihoodPID::GetPosterior(){
338 
339  return _posterior;
340 }
341 
342 double *LikelihoodPID::GetLikelihood(){
343 
344  return _likelihood;
345 }
346 
347 //track energy correction
348 double LikelihoodPID::getCorrEnergy(TLorentzVector pp, EVENT::Track* trk, EVENT::ClusterVec& cluvec){
349  int parttype=Classification(pp, trk, cluvec);
350  if(parttype==-1) parttype=2; //move to pion class
351 
352  return getCorrEnergy(pp, parttype);
353 }
354 
355 double LikelihoodPID::getCorrEnergy(TLorentzVector pp, int parttype){
356  if(parttype==-1) parttype=2; //move to pion class
357 
358  double tmpmass=0.0;
359  if(parttype==0) tmpmass=emass;
360  if(parttype==1) tmpmass=mmass;
361  if(parttype==2) tmpmass=pimass;
362  if(parttype==3) tmpmass=kmass;
363  if(parttype==4) tmpmass=pmass;
364 
365  return sqrt(pp.P()*pp.P()+tmpmass*tmpmass);
366 }
367 
368 //private
369 double LikelihoodPID::get_dEdxChi2(int parttype, TVector3 p, float hit, double dEdx){
370  //get parameters for chi2
371  double tmppar[5],tmpmass=0.0;
372  for(int i=0;i<5;i++) tmppar[i]=par[parttype][i];
373  //getmass
374  switch(parttype){
375  case 0:
376  tmpmass=emass;
377  break;
378  case 1:
379  tmpmass=mmass;
380  break;
381  case 2:
382  tmpmass=pimass;
383  break;
384  case 3:
385  tmpmass=kmass;
386  break;
387  case 4:
388  tmpmass=pmass;
389  break;
390  }
391 
392  //cal. polar angle
393  double trkcos=p.CosTheta();
394 
395  //get nomalized dEdx
396  double dEdx_Norm=get_Norm(dEdx);
397 
398  //get expected dE/dx
399  double ExpdEdx=BetheBloch(p.Mag(),tmpmass,tmppar);
400 
401  //cout << "check " << emass << " " << tmpmass << " " << trkcos << " " << dEdx_Norm << " " << ExpdEdx << endl;
402  double dEdx_Error = _dEdxerrfact * dEdx_Norm * hit/dEdx; //change 20151218
403 
404  //get chi2!!
405  double chi2=TMath::Power((dEdx_Norm-ExpdEdx)/dEdx_Error,2.0);
406  if(dEdx_Norm-ExpdEdx<0.0) chi2=-chi2; //get signed chi2
407 
408  return chi2;
409 }
410 
411 double LikelihoodPID::get_dEdxFactor(int parttype, TVector3 p, float hit, double dEdx){
412  //get parameters for chi2
413  double tmppar[5],tmpmass=0.0;
414  for(int i=0;i<5;i++) tmppar[i]=par[parttype][i];
415  //getmass
416  switch(parttype){
417  case 0:
418  tmpmass=emass;
419  break;
420  case 1:
421  tmpmass=mmass;
422  break;
423  case 2:
424  tmpmass=pimass;
425  break;
426  case 3:
427  tmpmass=kmass;
428  break;
429  case 4:
430  tmpmass=pmass;
431  break;
432  }
433 
434  //cal. polar angle
435  double trkcos=p.CosTheta();
436 
437  //get nomalized dEdx
438  double dEdx_Norm=get_Norm(dEdx);
439 
440  //cout << "check " << emass << " " << tmpmass << " " << trkcos << " " << dEdx_Norm << " " << ExpdEdx << endl;
441  double dEdx_Error = _dEdxerrfact * dEdx_Norm * hit/dEdx; //change 20151218
442 
443  //get likelihood factor
444  double factor=-0.5*TMath::Log(2.0*TMath::Pi()*dEdx_Error*dEdx_Error);
445 
446  return factor;
447 }
448 
449 //public
450 double LikelihoodPID::get_dEdxDist(int parttype){
451  //get parameters for dedxdist
452  double dedxdist=0.0;
453 
454  if(_dEdxDist[parttype]!=0.0) dedxdist = sqrt(fabs(_dEdxDist[parttype]))*_dEdxDist[parttype]/fabs(_dEdxDist[parttype]);
455  if(dedxdist!=dedxdist) dedxdist=0.0;
456  return dedxdist;
457 }
458 
459 double LikelihoodPID::get_Norm(double dedx){
460  return dedx/_dEdxnorm;
461 }
462 
463 double LikelihoodPID::BetheBloch(double x, double mass, double *pars){
464  double bg=x/mass;
465  double b=sqrt(bg*bg/(1.0+bg*bg));
466  //double g=bg/b;
467  double tmax=pars[2]*TMath::Power(bg,2.0); ///(1.0+pars[3]*g+pars[4]);
468 
469  return (0.5*pars[0]*TMath::Log(pars[1]*TMath::Power(bg,2.0)*tmax)-pars[3]*b*b-pars[4]*bg/2.0)/(b*b);
470 }
471 
472 int LikelihoodPID::Class_electron(TLorentzVector pp, EVENT::Track* trk, EVENT::ClusterVec& cluvec){
473  //Bayesian approach to classify electron!
474 
475  //get deposit energy
476  double ecal=0.0,hcal=0.0,mucal=0.0;
477  if(cluvec.size()!=0){
478  for(unsigned int i=0;i<cluvec.size();i++){
479  ecal+=cluvec[i]->getSubdetectorEnergies()[0];
480  hcal+=cluvec[i]->getSubdetectorEnergies()[1];
481  mucal+=cluvec[i]->getSubdetectorEnergies()[2];
482  }
483  }
484 
485  //get track shape parameters
486  if(cluvec.size()!=0) shapes=cluvec[0]->getShape();
487 
488  //get variables
489  double var[11]={}; //11
490  //for track variables
491  var[0]=(ecal+hcal)/pp.P();
492  if(ecal+hcal!=0.0) var[1]=ecal/(ecal+hcal);
493  if(shapes.size()!=0){
494  var[2]=shapes[0+4];
495  var[3]=shapes[5+4];
496  var[4]=fabs(shapes[3+4])/(shapes[6+4]);
497  var[5]=shapes[15+4]/(2.0*3.50);
498  }else var[2]=-1.0;
499 
500  var[6]=get_dEdxChi2(0,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
501  var[7]=get_dEdxChi2(1,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
502  var[8]=get_dEdxChi2(2,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
503  var[9]=get_dEdxChi2(3,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
504  var[10]=get_dEdxChi2(4,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
505 
506  var[6]=-0.5*fabs(var[6])+get_dEdxFactor(0,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
507  var[7]=-0.5*fabs(var[7])+get_dEdxFactor(1,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
508  var[8]=-0.5*fabs(var[8])+get_dEdxFactor(2,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
509  var[9]=-0.5*fabs(var[9])+get_dEdxFactor(3,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
510  var[10]=-0.5*fabs(var[10])+get_dEdxFactor(4,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
511 
512  //get likelihood for each class
513  double posterior[5]={0.0,0.0,0.0,0.0,0.0};
514  double risk[5]={0.0,0.0,0.0,0.0,0.0};
515  double okval[5]={0.0,0.0,0.0,0.0,0.0};
516  double total=0.0;
517  double priorprob[5]={0.0,0.0,0.0,0.0,0.0};
518  int valtype=0;
519  for(int j=0;j<7;j++){ //variables
520  //avoid very low momentum tracks (no energy deposit in the cal.)
521  if(var[0]==0.0 && j<=5) continue;
522 
523  //first, get likelihood
524  for(int i=0;i<5;i++){ //particle type
525  if(j==0) valtype=0;
526  if(j==1) valtype=1;
527  if(j==2) valtype=3;
528  if(j==3) valtype=4;
529  if(j==4) valtype=5;
530  if(j==5) valtype=6;
531  if(j==6) valtype=7;
532 
533  if(j<6){
534  okval[i]=getValue(i,valtype,var[j]); //likelihood
535  if(mucal==0.0 && i==1) okval[i]=getValue(5,valtype,var[j]); //likelihood
536  }else if(j==6)
537  okval[i]=std::max(std::exp(var[6+i]), 1.0e-300); //just use dE/dx likelihood
538  }
539 
540  //for basic variables flg
541  if(!_basicFlg && j<2) continue;
542  //for cluster shape flg
543  if(!_showerShapesFlg && (j>=2 && j<=5)) continue;
544  //for dEdx flg
545  if(!_dEdxFlg && j==6) continue;
546  //don't use some dEdxin the case of LikelihoodPID
547  if(_basicFlg && _showerShapesFlg && _dEdxFlg && !(j<=6)) continue;
548 
549  //to avoid strange value for dE/dx
550  //if(j>=6 && var[j]<=-50.0) continue;
551  if(j>=6 && var[j]!=var[j]) continue;
552 
553  //cal. prior probability
554  for(int i=0;i<5;i++) //particle type
555  priorprob[i]=prior[i]/
556  (prior[0]+prior[1]+prior[2]+prior[3]+prior[4]);
557 
558  //cal total probability for each particle type
559  total=0.0;
560  for(int i=0;i<5;i++) total+=priorprob[i]*okval[i];
561 
562  //cal. log posterior probability!
563  for(int i=0;i<5;i++){ //particle type
564  ///20150708 change it from posterior probability to likelihood
565  posterior[i]=TMath::Log(okval[i])+TMath::Log(priorprob[i])-TMath::Log(total);
566  _likelihood[i]+=TMath::Log(okval[i]);
567  }
568  //hadron type likelihood
569  _likelihood[5] += TMath::Log(_weights[2][valtype]*okval[2] + _weights[3][valtype]*okval[3] + _weights[4][valtype]*okval[4]);
570 
571  //bayesian updating
572  for(int i=0;i<5;i++) prior[i]=TMath::Exp(posterior[i]);
573 
574  }
575 
576  //cal. risk
577  //first. get penalty
578  // for(int i=0;i<5;i++){ //particle type
579  // for(int j=0;j<5;j++){
580  // penalty[i][j]=1.0;
581  // if(i==j) penalty[i][j]=1.0e-50;
582  // }
583  // }
584 
585  for(int i=0;i<5;i++){ //particle type
586  risk[i]=penalty[i][0]*TMath::Exp(posterior[0])+penalty[i][1]*TMath::Exp(posterior[1])
587  +penalty[i][2]*TMath::Exp(posterior[2])+penalty[i][3]*TMath::Exp(posterior[3])+penalty[i][4]*TMath::Exp(posterior[4]);
588  }
589 
590  //check the rule
591  int okflg=-1;
592  double tmppp=-1.0e+100;
593  if(_usebayes==0){
594  //---------- here is likelihood based determination rule -----------
595  for(int i=0;i<5;i++){
596  tmppp=TMath::Max(_likelihood[i], tmppp);
597  //if(posterior[0]<posterior[i]){
598  // okflg=false;
599  // break;
600  //}
601  }
602 
603  for(int i=0;i<5;i++){
604  //if(fabs(tmppp-posterior[i])<1.0e-6) okflg=i;
605  if(fabs(tmppp-_likelihood[i])<1.0e-6) okflg=i;
606  }
607  //------------------------------------------------------------------
608  }else{
609  //---------- here is risk based determination rule -----------
610  tmppp=1.0e+100;
611  for(int i=0;i<5;i++){
612  tmppp=TMath::Min(risk[i], tmppp);
613  }
614 
615  for(int i=0;i<5;i++){
616  if(fabs(tmppp-risk[i])<1.0e-8) okflg=i;
617  }
618 
619  //threshold is set.
620  if(okflg==0 && posterior[okflg]<TMath::Log(threshold[okflg])) okflg=-1;
621  //------------------------------------------------------------------
622  }
623 
624  //save posterior
625  _posterior[0]=TMath::Exp(posterior[0]);
626  _posterior[1]=TMath::Exp(posterior[1]);
627  _posterior[2]=TMath::Exp(posterior[2]);
628  _posterior[3]=TMath::Exp(posterior[3]);
629  _posterior[4]=TMath::Exp(posterior[4]);
630  //hadron type
631  _posterior[5]=_posterior[2]+_posterior[3]+_posterior[4];
632 
633  return okflg;
634 }
635 
636 int LikelihoodPID::Class_muon(TLorentzVector pp, EVENT::Track* trk, EVENT::ClusterVec& cluvec){
637  //Bayesian approach to classify muon!
638 
639  //get deposit energy
640  double ecal=0.0,hcal=0.0,mucal=0.0;
641  if(cluvec.size()!=0){
642  for(unsigned int i=0;i<cluvec.size();i++){
643  ecal+=cluvec[i]->getSubdetectorEnergies()[0];
644  hcal+=cluvec[i]->getSubdetectorEnergies()[1];
645  mucal+=cluvec[i]->getSubdetectorEnergies()[2];
646  }
647  }
648 
649  //get track shape parameters
650  if(cluvec.size()!=0) shapes=cluvec[0]->getShape();
651 
652  //get variables
653  double var[17]={}; //17 variales used(is it OK?)
654  var[0]=(ecal+hcal)/pp.P();
655  if(ecal+hcal!=0.0) var[1]=ecal/(ecal+hcal);
656  var[2]=mucal;
657 
658  if(shapes.size()!=0){
659  var[3]=shapes[0+4];
660  var[4]=shapes[5+4];
661  var[5]=fabs(shapes[3+4])/(shapes[6+4]);
662  var[6]=shapes[15+4]/(2.0*3.50);
663 
664  //for low momentum mu/pi separation
665  var[14]= shapes[17+4];
666  var[15]= shapes[18+4];
667  var[16]= shapes[19+4];
668  }else var[3]=-1.0;
669 
670  var[7]=get_dEdxChi2(0,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
671  var[8]=get_dEdxChi2(1,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
672  var[9]=get_dEdxChi2(2,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
673  var[10]=get_dEdxChi2(3,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
674  var[11]=get_dEdxChi2(4,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
675 
676  var[7]=-0.5*fabs(var[7])+get_dEdxFactor(0,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
677  var[8]=-0.5*fabs(var[8])+get_dEdxFactor(1,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
678  var[9]=-0.5*fabs(var[9])+get_dEdxFactor(2,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
679  var[10]=-0.5*fabs(var[10])+get_dEdxFactor(3,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
680  var[11]=-0.5*fabs(var[11])+get_dEdxFactor(4,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
681 
682  var[12] = _delpos[0];
683  var[13] = _delpos[2];
684 
685  //get likelihood for each class
686  double posterior[5]={0.0,0.0,0.0,0.0,0.0};
687  double risk[5]={0.0,0.0,0.0,0.0,0.0};
688  int valtype=0;
689  double okval[5]={0.0,0.0,0.0,0.0,0.0};
690  double total=0.0; //[5]={0.0,0.0,0.0,0.0,0.0};
691  double priorprob[5]={0.0,0.0,0.0,0.0,0.0};
692  for(int j=0;j<13;j++){ //variables
693 
694  if(var[0]>0.0){
695  //avoid pion misID when energy deposit to mucal is zero
696  if(var[2]==0.0 && j==2) continue; //(j==2 || j==7 || j>=9)) continue;
697  //don't use when ep>0.0
698  //else if(var[2]>0.0 && j>=9) continue;
699  }else{
700  //avoid very low momentum tracks (no energy deposit in the cal.)
701  if(j<=6) continue;
702  if(j>7) continue;
703  }
704 
705  //first, get likelihood
706  for(int i=0;i<5;i++){ //particle type(electron isn't be checked)
707  if(j==0) valtype=0;
708  if(j==1) valtype=1;
709  if(j==2) valtype=2;
710  if(j==3) valtype=3;
711  if(j==4) valtype=4;
712  if(j==5) valtype=5;
713  if(j==6) valtype=6;
714  if(j==7) valtype=7;
715  if(j==8) valtype=11;
716  if(j==9) valtype=12;
717  if(j==10) valtype=13;
718  if(j==11) valtype=14;
719  if(j==12) valtype=15;
720 
721  if(j<7){
722  okval[i]=getValue(i,valtype,var[j]); //likelihood
723  if(var[2]==0.0 && i==1) okval[i]=getValue(5,valtype,var[j]); //likelihood
724  }else if(j==7){
725  okval[i]=std::max(std::exp(var[7+i]), 1.0e-300); //just use dE/dx likelihood
726  }else if(j==8 || j==9){
727  okval[i]=getValue(i,valtype,var[j+4]); //likelihood
728  if(var[2]==0.0 && i==1) okval[i]=getValue(5,valtype,var[j+4]); //likelihood
729  }else if(j>=10 && pp.P()>0.3 && pp.P()<6.0){ //for low momentum mu/pi
730  okval[i]=getValue(i,valtype,var[j+4]); //likelihood
731  if(var[2]==0.0 && i==1) okval[i]=getValue(5,valtype,var[j+4]); //likelihood
732  }
733  }
734 
735  //for basic variables flg
736  if(!_basicFlg && j<=2) continue;
737  //for cluster shape flg
738  if(!_showerShapesFlg && ((j>2 && j<=6) || j>=8)) continue;
739  //for dEdx flg
740  if(!_dEdxFlg && j==7) continue;
741  //don't use some dEdxin the case of LikelihoodPID
742  if(_basicFlg && _showerShapesFlg && _dEdxFlg && !(j<=9 || (j>=10 && pp.P()>0.3 && pp.P()<6.0))) continue;
743  if(mucal!=0.0 && j==10) continue;
744 
745  //to avoid strange value for dE/dx
746  //if(j>=7 && j<=11 && var[j]<=-50.0) continue;
747  if(j>=7 && j<=11 && var[j]!=var[j]) continue;
748 
749  //cal. prior probability
750  for(int i=0;i<5;i++) //particle type
751  priorprob[i]=prior[i]/
752  (prior[0]+prior[1]+prior[2]+prior[3]+prior[4]);
753 
754  //cal total probability for each particle type
755  total=0.0;
756  for(int i=0;i<5;i++) total+=priorprob[i]*okval[i];
757 
758  //cal. log likelihood first
759  for(int i=0;i<5;i++){ //particle type
760  //20150708 change it from posterior probability to simple likelihood
761  posterior[i]=TMath::Log(okval[i])+TMath::Log(priorprob[i])-TMath::Log(total);
762  _likelihood[i]+=TMath::Log(okval[i]);
763  }
764  //hadron type likelihood
765  _likelihood[5] += TMath::Log(_weights[2][valtype]*okval[2] + _weights[3][valtype]*okval[3] + _weights[4][valtype]*okval[4]);
766 
767  //bayesian updating
768  for(int i=0;i<5;i++) prior[i]=TMath::Exp(posterior[i]);
769  }
770 
771  //cal. risk
772  //first. get penalty
773  // for(int i=0;i<5;i++){ //particle type
774  // for(int j=0;j<5;j++){
775  // penalty[i][j]=1.0;
776  // if(i==j) penalty[i][j]=1.0e-50;
777  // }
778  // }
779  // penalty[1][2] = 9.0; //TMath::Min(1.0, 9.0*pp.P(); Todo: momentum depepndence
780 
781  for(int i=0;i<5;i++){ //particle type
782  risk[i]=penalty[i][0]*TMath::Exp(posterior[0])+penalty[i][1]*TMath::Exp(posterior[1])
783  +penalty[i][2]*TMath::Exp(posterior[2])+penalty[i][3]*TMath::Exp(posterior[3])+penalty[i][4]*TMath::Exp(posterior[4]);
784  }
785 
786  //check the rule
787  int okflg=-1;
788  double tmppp=-1.0e+100;
789  if(_usebayes==0){
790  //---------- here is likelihood based determination rule -----------
791  for(int i=0;i<5;i++){
792  tmppp=TMath::Max(_likelihood[i], tmppp);
793  //if(posterior[0]<posterior[i]){
794  // okflg=false;
795  // break;
796  //}
797  }
798 
799  for(int i=0;i<5;i++){
800  //if(fabs(tmppp-posterior[i])<1.0e-6) okflg=i;
801  if(fabs(tmppp-_likelihood[i])<1.0e-6) okflg=i;
802  }
803  //------------------------------------------------------------------
804  }else{
805  //---------- here is risk based determination rule -----------
806  tmppp=1.0e+100;
807  for(int i=0;i<5;i++){
808  tmppp=TMath::Min(risk[i], tmppp);
809  }
810 
811  for(int i=0;i<5;i++){
812  if(fabs(tmppp-risk[i])<1.0e-8) okflg=i;
813  }
814 
815  if(okflg==1 && posterior[okflg]<TMath::Log(threshold[okflg])) okflg=-1; //threshold is set.
816  //------------------------------------------------------------------
817  }
818 
819  //save posterior
820  _posterior[0]=TMath::Exp(posterior[0]);
821  _posterior[1]=TMath::Exp(posterior[1]);
822  _posterior[2]=TMath::Exp(posterior[2]);
823  _posterior[3]=TMath::Exp(posterior[3]);
824  _posterior[4]=TMath::Exp(posterior[4]);
825  //hadron type
826  _posterior[5]=_posterior[2]+_posterior[3]+_posterior[4];
827 
828  return okflg;
829 }
830 
831 int LikelihoodPID::Class_hadron(TLorentzVector pp, EVENT::Track* trk, EVENT::ClusterVec& cluvec){
832  //Bayesian approach to classify hadrons!
833 
834  //get deposit energy
835  double ecal=0.0,hcal=0.0,mucal=0.0;
836  if(cluvec.size()!=0){
837  for(unsigned int i=0;i<cluvec.size();i++){
838  ecal+=cluvec[i]->getSubdetectorEnergies()[0];
839  hcal+=cluvec[i]->getSubdetectorEnergies()[1];
840  mucal+=cluvec[i]->getSubdetectorEnergies()[2];
841  }
842  }
843 
844  //get track shape parameters
845  if(cluvec.size()!=0) shapes=cluvec[0]->getShape();
846 
847  //get variables
848  double var[10]; //10 variales used(is it OK?)
849  var[0]=(ecal+hcal)/pp.P();
850  if(ecal+hcal!=0.0)var[1]=ecal/(ecal+hcal);
851  if(shapes.size()!=0){
852  var[2]=shapes[0+4];
853  var[3]=fabs(shapes[3+4])/(shapes[6+4]);
854  var[4]=shapes[15+4]/(2.0*3.50);
855  }else var[2]=-1.0;
856 
857  var[5]=get_dEdxChi2(0,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
858  var[6]=get_dEdxChi2(1,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
859  var[7]=get_dEdxChi2(2,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
860  var[8]=get_dEdxChi2(3,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
861  var[9]=get_dEdxChi2(4,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //chi2
862 
863  var[5]=-0.5*fabs(var[5])+get_dEdxFactor(0,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
864  var[6]=-0.5*fabs(var[6])+get_dEdxFactor(1,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
865  var[7]=-0.5*fabs(var[7])+get_dEdxFactor(2,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
866  var[8]=-0.5*fabs(var[8])+get_dEdxFactor(3,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
867  var[9]=-0.5*fabs(var[9])+get_dEdxFactor(4,pp.Vect(),trk->getdEdxError(),trk->getdEdx()); //log likelihood
868 
869  //get likelihood for each class
870  double posterior[5]={TMath::Log(0.20),TMath::Log(0.20),TMath::Log(0.20),TMath::Log(0.20),TMath::Log(0.20)};
871  double risk[5]={0.0,0.0,0.0,0.0,0.0};
872  int valtype=0;
873  double okval[5]={0.0,0.0,0.0,0.0,0.0};
874  double total[5]={0.0,0.0,0.0,0.0,0.0};
875  double priorprob[5]={0.0,0.0,0.0,0.0,0.0};
876  for(int j=0;j<6;j++){ //variables
877  //avoid very low momentum tracks (no energy deposit in the cal.)
878  if(var[0]==0.0 && j<=4) continue;
879 
880  //first, get likelihood
881  for(int i=0;i<5;i++){ //particle type(electron&muon aren't be checked)
882  if(j==0) valtype=0;
883  if(j==1) valtype=1;
884  if(j==2) valtype=3;
885  if(j==3) valtype=5;
886  if(j==4) valtype=6;
887  if(j==5) valtype=7;
888  if(j==6) valtype=8;
889  if(j==7) valtype=9;
890  if(j==8) valtype=10;
891  if(j==9) valtype=11;
892 
893  if(j<5){
894  okval[i]=getValue(i,valtype,var[j]); //likelihood
895  }else if(j>=5)
896  okval[i]=std::max(std::exp(var[5+i]), 1.0e-300); //just use dE/dx likelihood
897  }
898 
899  //for basic variables flg
900  if(!_basicFlg && j<2) continue;
901  //for cluster shape flg
902  if(!_showerShapesFlg && j>=2 && j<=4) continue;
903  //for dEdx flg
904  if(!_dEdxFlg && j>=5) continue;
905  //don't use some dEdxin the case of LikelihoodPID
906  if(_basicFlg && _showerShapesFlg && _dEdxFlg && !(j<=5)) continue;
907 
908  //to avoid strange value for dE/dx
909  //if(j>=5 && var[j]<=-50.0) continue;
910  if(j>=5 && var[5]!=var[5]) continue;
911 
912  //cal. prior probability
913  for(int i=0;i<5;i++) //particle type
914  priorprob[i]=prior[i]/
915  (prior[0]+prior[1]+prior[2]+prior[3]+prior[4]);
916 
917  //cal total probability
918  for(int i=0;i<5;i++){ //particle type
919  total[i]=priorprob[0]*okval[0]+priorprob[1]*okval[1]+priorprob[2]*okval[2]+priorprob[3]*okval[3]+priorprob[4]*okval[4];
920  if(total[i]==0.0) total[i]=1.0e-100;
921  }
922 
923  //cal. log posterior probability!
924  for(int i=0;i<5;i++){ //particle type
925  //20150708 change it from posterior probabity to simple likelihood
926  posterior[i]=TMath::Log(okval[i])+TMath::Log(priorprob[i])-TMath::Log(total[i]);
927  _likelihood[i]+=TMath::Log(okval[i]);
928  }
929  //hadron type likelihood
930  _likelihood[5] += TMath::Log(_weights[2][valtype]*okval[2] + _weights[3][valtype]*okval[3] + _weights[4][valtype]*okval[4]);
931 
932  //bayesian updating
933  for(int i=0;i<5;i++) prior[i]=TMath::Exp(posterior[i]);
934  }
935 
936  //cal. risk
937  //normalize
938  // for(int i=0;i<5;i++){
939  // for(int j=0;j<5;j++){
940  // penalty[i][j]=1.0;
941  // if(i==j) penalty[i][j]=1.0e-50;
942  // }
943  // }
944  //first. get penalty
945  // for(int i=2;i<5;i++){ //particle type
946  // for(int j=2;j<5;j++){
947  // penalty[i][j]=getPenalty(i,j,pp.P());
948  // }
949  // double tmppenalty=penalty[i][i];
950  // for(int j=2;j<5;j++){
951  // penalty[i][j]=(fabs(penalty[i][j]-tmppenalty)+tmppenalty)/tmppenalty;
952  // penalty[i][j]=penalty[i][j]-1.0;
953 
954  // }
955  // }
956 
957  //weight of the cost
958  //reserve default value
959  double res43=penalty[4][3];
960  double res32=penalty[3][2];
961  if(_usebayes==2){
962  if(pp.P()>1.0){
963  penalty[4][3] = max(1.0, penalty[4][3]*(1.0-TMath::Exp(-0.772*pp.P()))); //10.0
964  //penalty[3][2] = max(1.0, penalty[3][2]*TMath::Exp(0.0268*pp.P())); //30.0
965  }else{
966  penalty[4][3] = 1.0;
967  penalty[3][2] = 1.0;
968  }
969 
970  if(fabs(pp.CosTheta())<1.0){
971  penalty[4][3] = max(1.0, penalty[4][3]*(0.439*pp.CosTheta()*pp.CosTheta()-0.940*fabs(pp.CosTheta())+0.495));
972  penalty[3][2] = max(1.0, penalty[3][2]*TMath::Exp(-1.46*fabs(pp.CosTheta()))); //30.0
973  }
974 
975  }
976 
977  //penalty[4][2] = 5.0;
978 
979  for(int i=0;i<5;i++){ //particle type
980  risk[i]=penalty[i][0]*TMath::Exp(posterior[0])+penalty[i][1]*TMath::Exp(posterior[1])
981  +penalty[i][2]*TMath::Exp(posterior[2])+penalty[i][3]*TMath::Exp(posterior[3])+penalty[i][4]*TMath::Exp(posterior[4]);
982  }
983 
984  penalty[4][3]=res43;
985  penalty[3][2]=res32;
986 
987  //check the rule
988  int okflg=-1;
989  double tmppp=-1.0e+100;
990  if(_usebayes==0){
991  //---------- here is likelihood based determination rule -----------
992  for(int i=0;i<5;i++){
993  tmppp=TMath::Max(_likelihood[i], tmppp);
994  //if(posterior[0]<posterior[i]){
995  // okflg=false;
996  // break;
997  //}
998  }
999 
1000  for(int i=0;i<5;i++){
1001  //if(fabs(tmppp-posterior[i])<1.0e-6) okflg=i;
1002  if(fabs(tmppp-_likelihood[i])<1.0e-6) okflg=i;
1003  }
1004 
1005  //move to pion when very bad
1006  //if(TMath::Exp(tmppp)<0.21) okflg=2;
1007  if(tmppp==0.0) okflg=2;
1008  //------------------------------------------------------------------
1009  }else{
1010  //---------- here is risk based determination rule -----------
1011  tmppp=1.0e+100;
1012  for(int i=0;i<5;i++){
1013  tmppp=TMath::Min(risk[i], tmppp);
1014  }
1015 
1016  for(int i=0;i<5;i++){
1017  if(fabs(tmppp-risk[i])<1.0e-6) okflg=i;
1018  }
1019  if(okflg<2) okflg=2;
1020  if(posterior[okflg]<=TMath::Log(0.21)) okflg=TMath::Max(2,okflg); //do not go to leptons
1021 
1022  //check penalty -1 means undefined
1023  if(okflg>=2 && posterior[okflg]<TMath::Log(threshold[okflg])) okflg=-1; //threshold is set.
1024  //------------------------------------------------------------------
1025  }
1026 
1027  //save posterior
1028  _posterior[0]=TMath::Exp(posterior[0]);
1029  _posterior[1]=TMath::Exp(posterior[1]);
1030  _posterior[2]=TMath::Exp(posterior[2]);
1031  _posterior[3]=TMath::Exp(posterior[3]);
1032  _posterior[4]=TMath::Exp(posterior[4]);
1033  //hadron type
1034  _posterior[5]=_posterior[2]+_posterior[3]+_posterior[4];
1035 
1036  return okflg;
1037 }
1038 
1039 double LikelihoodPID::getValue(int type, int valtype, double value){
1040 
1041  /*int nbins=pdf[type][valtype]->GetNbinsX();
1042  double interval=pdf[type][valtype]->GetBinWidth(1);
1043  double minval=pdf[type][valtype]->GetBinLowEdge(1);
1044 
1045  double val=1.0e-30;
1046  int bin=0;
1047  for(int i=0;i<nbins;i++){
1048  if(valtype!=1 && valtype<7){
1049  if(value<minval){
1050  bin=0;
1051  break;
1052  }
1053  if(value>=minval+i*interval && value<minval+(i+1)*interval){
1054  bin=i+1;
1055  break;
1056  }
1057  if(value>=minval+nbins*interval){
1058  bin=nbins+1;
1059  break;
1060  }
1061  }else if(valtype>=7){ //for hadem distribution!(because 1.0 is the max value!)
1062  if(value<=minval){
1063  bin=0;
1064  break;
1065  }
1066  if(value>minval+i*interval && value<=minval+(i+1)*interval){
1067  bin=i+1;
1068  break;
1069  }
1070  if(value>minval+nbins*interval){
1071  bin=nbins+1;
1072  break;
1073  }
1074  }else if(valtype==1){ //ehad very strange graph!
1075  if(value<=minval){
1076  bin=1;
1077  break;
1078  }
1079  if(value>minval+i*interval && value<=minval+(i+1)*interval){
1080  bin=i+1;
1081  break;
1082  }
1083  if(value>minval+nbins*interval){
1084  bin=nbins+1;
1085  break;
1086  }
1087  }
1088  }*/
1089 
1090 
1091  double val=1.0e-100;
1092  int bin=0;
1093 
1094  if(pdf[type][valtype]!=NULL){
1095  bin=pdf[type][valtype]->GetXaxis()->FindBin(value);
1096 
1097  //get probability
1098  val=pdf[type][valtype]->GetBinContent(bin);
1099  if(val==0.0) val= 1.0e-100;
1100  }else
1101  val = 1.0;
1102 
1103  return val;
1104 }
1105 
1106 double LikelihoodPID::getPenalty(int ptype, int hypothesis, double p){
1107  double par[3]={0.0,0.0,0.0};
1108  //set parameters
1109  switch(ptype){
1110  case 0: //electron
1111  if(hypothesis==0){
1112  par[0]=0.0;
1113  par[1]=1.0;
1114  par[2]=1.00178;
1115 
1116  }else if(hypothesis==1){
1117  par[0]=0.101945;
1118  par[1]=2.18719;
1119  par[2]=0.990000;
1120 
1121  }else if(hypothesis==2){
1122  par[0]=0.104521;
1123  par[1]=1.98490;
1124  par[2]=0.990000;
1125 
1126  }else if(hypothesis==3){
1127  par[0]=0.234299;
1128  par[1]=0.276835;
1129  par[2]=0.965973;
1130 
1131  }else{
1132  par[0]=0.414613;
1133  par[1]=-0.132530;
1134  par[2]=0.949679;
1135 
1136  }
1137 
1138  break;
1139  case 1: //muon
1140  if(hypothesis==0){
1141  par[0]=-0.00512190;
1142  par[1]=-0.211835;
1143  par[2]=1.00024;
1144 
1145  }else if(hypothesis==1){
1146  par[0]=0.0;
1147  par[1]=1.00;
1148  par[2]=0.999684;
1149 
1150  }else if(hypothesis==2){
1151  par[0]=0.00833283;
1152  par[1]=9.99995;
1153  par[2]=0.999027;
1154 
1155  }else if(hypothesis==3){
1156  par[0]=0.0964021;
1157  par[1]=-0.214469;
1158  par[2]=0.989688;
1159 
1160  }else{
1161  par[0]=0.318674;
1162  par[1]=-0.197755;
1163  par[2]=0.968436;
1164  }
1165 
1166  break;
1167  case 2: //pion
1168  if(hypothesis==0){
1169  par[0]=-0.0123577;
1170  par[1]=-0.141521;
1171  par[2]=1.00273;
1172 
1173  }else if(hypothesis==1){
1174  par[0]=-0.00558462;
1175  par[1]=-0.136941;
1176  par[2]=1.00135;
1177 
1178  }else if(hypothesis==2){
1179  par[0]=0.0;
1180  par[1]=1.0;
1181  par[2]=1.00001;
1182 
1183  }else if(hypothesis==3){
1184  par[0]=0.122083;
1185  par[1]=-0.1333923;
1186  par[2]=0.976863;
1187 
1188  }else{
1189  par[0]=0.401111;
1190  par[1]=-0.116807;
1191  par[2]=0.930906;
1192  }
1193 
1194  break;
1195  case 3: //kaon
1196  if(hypothesis==0){
1197  par[0]=-0.102300;
1198  par[1]=-0.139570;
1199  par[2]=1.01362;
1200 
1201  }else if(hypothesis==1){
1202  par[0]=-0.0973257;
1203  par[1]=-0.138932;
1204  par[2]=1.01293;
1205 
1206  }else if(hypothesis==2){
1207  par[0]=-0.0936450;
1208  par[1]=-0.138469;
1209  par[2]=1.01242;
1210 
1211  }else if(hypothesis==3){
1212  par[0]=0.0;
1213  par[1]=1.0;
1214  par[2]=0.999865;
1215 
1216  }else{
1217  par[0]=0.223317;
1218  par[1]=-0.101273;
1219  par[2]=0.973484;
1220  }
1221 
1222  break;
1223  case 4: //proton
1224  if(hypothesis==0){
1225  par[0]=-0.260150;
1226  par[1]=-0.0990612;
1227  par[2]=1.02854;
1228 
1229  }else if(hypothesis==1){
1230  par[0]=-0.256503;
1231  par[1]=-0.0976728;
1232  par[2]=1.02811;
1233 
1234  }else if(hypothesis==2){
1235  par[0]=-0.253788;
1236  par[1]=-0.0966732;
1237  par[2]=1.02779;
1238 
1239  }else if(hypothesis==3){
1240  par[0]=-0.183031;
1241  par[1]=-0.0742123;
1242  par[2]=1.01965;
1243 
1244  }else{
1245  par[0]=0.0;
1246  par[1]=1.0;
1247  par[2]=0.999791;
1248  }
1249 
1250  break;
1251  }
1252 
1253  return par[0]/sqrt(p*p+par[1])+par[2];
1254 }
1255 
1256 void LikelihoodPID::CalculateDeltaPosition(float charge, TVector3 p, const float* calpos){
1257  //calculate some variables for preparation
1258  //transverse momentum
1259  double prphi=sqrt(pow(p[0],2)+pow(p[1],2));
1260  //momentum
1261  double pp=p.Mag();
1262  //radius in r-phi plane
1263  double radius=prphi/(0.3*_bfield);
1264  //tangent lambda
1265  double tanlam=p[2]/prphi;
1266  //change to unit vector
1267  p[0]=p[0]/prphi;
1268  p[1]=p[1]/prphi;
1269  p[2]=p[2]/pp;
1270 
1271  //chenge position vector to meter
1272  float calcalpos[3];
1273  calcalpos[0]=calpos[0]/1000.0;
1274  calcalpos[1]=calpos[1]/1000.0;
1275  calcalpos[2]=calpos[2]/1000.0;
1276 
1277 
1278  //radius at the tracker end
1279  double rradius=sqrt(pow(calcalpos[0],2)+pow(calcalpos[1],2));
1280 
1281  //cout << "check val. " << prphi << " " << radius << " " << tanlam
1282  // << " " << rradius << endl;
1283 
1284  //cal. the position of the center of track circle
1285  TVector3 cc;
1286  cc[0]=-charge*p[1]*prphi/(0.3*_bfield);
1287  cc[1]=charge*p[0]*prphi/(0.3*_bfield);
1288  cc[2]=radius*tanlam;
1289 
1290  //cal. sign and cosine 2theta
1291  double sintheta=charge*rradius/(2*radius);
1292  double costheta=sqrt(1-sintheta*sintheta);
1293  double sin2theta=2*sintheta*costheta;
1294  double cos2theta=costheta*costheta-sintheta*sintheta;
1295 
1296  TVector3 calpos2;
1297  calpos2[2]=rradius*tanlam;
1298  calpos2[0]=cc[0]*(1-cos2theta)+cc[1]*sin2theta;
1299  calpos2[1]=-cc[0]*sin2theta+cc[1]*(1-cos2theta);
1300 
1301  //cal. difference of the position
1302  //float delpos[3];
1303  _delpos[0]=charge*fabs(calcalpos[0]-calpos2[0]);
1304  _delpos[1]=charge*fabs(calcalpos[1]-calpos2[1]);
1305  _delpos[2]=fabs(calcalpos[2]-calpos2[2]);
1306 
1307 
1308  //cout << "calpos: " << calcalpos[0] << " " << calcalpos[1] << " " << calcalpos[2] << endl;
1309  //cout << "calpos2: " << calpos2[0] << " " << calpos2[1] << " " << calpos2[2] << endl;
1310  //cout << "result: " << delpos[0] << " " << delpos[1] << " " << delpos[2] << endl;
1311  return;
1312 }
1313 
LikelihoodPID::Class_hadron
int Class_hadron(TLorentzVector pp, EVENT::Track *trk, EVENT::ClusterVec &cluvec)
Definition: LikelihoodPID.cc:831
LikelihoodPID::get_dEdxDist
double get_dEdxDist(int parttype)
Definition: LikelihoodPID.cc:450
LikelihoodPID::CalculateDeltaPosition
void CalculateDeltaPosition(float charge, TVector3 p, const float *calpos)
Definition: LikelihoodPID.cc:1256
LikelihoodPID::Class_muon
int Class_muon(TLorentzVector pp, EVENT::Track *trk, EVENT::ClusterVec &cluvec)
Definition: LikelihoodPID.cc:636
LikelihoodPID::Class_electron
int Class_electron(TLorentzVector pp, EVENT::Track *trk, EVENT::ClusterVec &cluvec)
Definition: LikelihoodPID.cc:472
LikelihoodPID::GetLikelihood
double * GetLikelihood()
Definition: LikelihoodPID.cc:342
LikelihoodPID::get_dEdxChi2
double get_dEdxChi2(int parttype, TVector3 p, float hit, double dEdx)
Definition: LikelihoodPID.cc:369
LikelihoodPID::get_dEdxFactor
double get_dEdxFactor(int parttype, TVector3 p, float hit, double dEdx)
Definition: LikelihoodPID.cc:411
LikelihoodPID::GetPosterior
double * GetPosterior()
Definition: LikelihoodPID.cc:337
sistrip::s
static const float s
Definition: PhysicalConstants.h:50
LikelihoodPID::getValue
double getValue(int type, int valtype, double value)
Definition: LikelihoodPID.cc:1039
itos
std::string itos(int i)
Definition: LikelihoodPID.cc:34
sistrip::e
static const float e
Definition: PhysicalConstants.h:26
LikelihoodPID::get_Norm
double get_Norm(double dedx)
Definition: LikelihoodPID.cc:459
LikelihoodPID::getCorrEnergy
double getCorrEnergy(TLorentzVector pp, EVENT::Track *trk, EVENT::ClusterVec &cluvec)
Definition: LikelihoodPID.cc:348
LikelihoodPID::Classification
int Classification(TLorentzVector pp, EVENT::Track *trk, EVENT::ClusterVec &cluvec)
Definition: LikelihoodPID.cc:277
LikelihoodPID::getPenalty
double getPenalty(int ptype, int hypothesis, double p)
Definition: LikelihoodPID.cc:1106
LikelihoodPID::BetheBloch
double BetheBloch(double x, double mass, double *pars)
Definition: LikelihoodPID.cc:463