CbmRoot
CbmKresTrainAnn.cxx
Go to the documentation of this file.
1 
18 #include "CbmKresTrainAnn.h"
19 
20 #include <boost/assign/list_of.hpp>
21 
22 #include <cmath>
23 #include <iostream>
24 #include <string>
25 #include <vector>
26 
27 #include "CbmDrawHist.h"
28 #include "TCanvas.h"
29 #include "TH1D.h"
30 #include "TH2D.h"
31 #include "TMath.h"
32 #include "TMultiLayerPerceptron.h"
33 #include "TSystem.h"
34 #include "TTree.h"
35 
36 
37 using boost::assign::list_of;
38 using namespace std;
39 
40 CbmKresTrainAnn::CbmKresTrainAnn()
41  : fMaxNofTrainSamples(1500)
42  , fAnnCut(0.6)
43  , fNofWrongLikeCorrect(0)
44  , fNofCorrectLikeWrong(0)
45  , IM_correct()
46  , OA_correct()
47  , Angle_correct()
48  , Z_correct()
49  , Mom1_correct()
50  , Mom2_correct()
51  , IM_wrong()
52  , OA_wrong()
53  , Angle_wrong()
54  , Z_wrong()
55  , Mom1_wrong()
56  , Mom2_wrong()
57  , fHists()
58  , fhAnnOutput_correct(nullptr)
59  , fhAnnOutput_wrong(nullptr)
60  , fhCumProb_correct(nullptr)
61  , fhCumProb_wrong(nullptr) {}
62 
63 CbmKresTrainAnn::~CbmKresTrainAnn() {}
64 
65 void CbmKresTrainAnn::Init() { InitHistograms(); }
66 
67 void CbmKresTrainAnn::Exec(int event,
68  int IdForANN,
69  double InvariantMass,
70  double OpeningAngle,
71  double PlaneAngle_last,
72  double ZPos,
73  TVector3 Momentum1,
74  TVector3 Momentum2) {
75  double p1 =
76  TMath::Sqrt(Momentum1.X() * Momentum1.X() + Momentum1.Y() * Momentum1.Y()
77  + Momentum1.Z() * Momentum1.Z());
78  double p2 =
79  TMath::Sqrt(Momentum2.X() * Momentum2.X() + Momentum2.Y() * Momentum2.Y()
80  + Momentum2.Z() * Momentum2.Z());
81  if (IdForANN == 1) {
82  //if (IM_correct.size() < fMaxNofTrainSamples){
83  IM_correct.push_back(InvariantMass);
84  OA_correct.push_back(OpeningAngle);
85  Angle_correct.push_back(PlaneAngle_last);
86  Z_correct.push_back(ZPos);
87  Mom1_correct.push_back(p1);
88  Mom2_correct.push_back(p2);
89  //}
90  } else {
91  //if (IM_wrong.size() < fMaxNofTrainSamples){
92  IM_wrong.push_back(InvariantMass);
93  OA_wrong.push_back(OpeningAngle);
94  Angle_wrong.push_back(PlaneAngle_last);
95  Z_wrong.push_back(ZPos);
96  Mom1_wrong.push_back(p1);
97  Mom2_wrong.push_back(p2);
98  //}
99  }
100 
101  if (IM_correct.size() % 100 == 0 && IdForANN == 1)
102  cout << "correct = " << IM_correct.size()
103  << "; wrong = " << IM_wrong.size() << endl;
104 
105  //if (IM_correct.size() >= fMaxNofTrainSamples && IM_wrong.size() >= fMaxNofTrainSamples) {
106  if (event == 2000 && IM_correct.size() >= fMaxNofTrainSamples) {
107  TrainAndTestAnn();
108  Draw();
109 
110  IM_correct.clear();
111  OA_correct.clear();
112  Angle_correct.clear();
113  Z_correct.clear();
114  Mom1_correct.clear();
115  Mom2_correct.clear();
116  IM_wrong.clear();
117  OA_wrong.clear();
118  Angle_wrong.clear();
119  Z_wrong.clear();
120  Mom1_wrong.clear();
121  Mom2_wrong.clear();
122  }
123 }
124 
125 void CbmKresTrainAnn::TrainAndTestAnn() {
126  //cout << "Do TrainAndTestAnn" << endl;
127  TTree* simu = new TTree("MonteCarlo", "MontecarloData");
128  Double_t x[6];
129  Double_t xOut;
130 
131  simu->Branch("x0", &x[0], "x0/D");
132  simu->Branch("x1", &x[1], "x1/D");
133  simu->Branch("x2", &x[2], "x2/D");
134  simu->Branch("x3", &x[3], "x3/D");
135  simu->Branch("x4", &x[4], "x4/D");
136  simu->Branch("x5", &x[5], "x5/D");
137  simu->Branch("xOut", &xOut, "xOut/D");
138 
139  for (size_t i = 0; i < IM_correct.size(); i++) {
140  x[0] = IM_correct[i] / 0.1;
141  x[1] = OA_correct[i] / 30;
142  x[2] = Angle_correct[i] / 30;
143  x[3] = Z_correct[i] / 100;
144  x[4] = Mom1_correct[i] / 5;
145  x[5] = Mom2_correct[i] / 5;
146 
147  if (x[0] > 1.0) x[0] = 1.0;
148  if (x[1] > 1.0) x[1] = 1.0;
149  if (x[2] > 1.0) x[2] = 1.0;
150  if (x[3] > 1.0) x[3] = 1.0;
151  if (x[4] > 1.0) x[4] = 1.0;
152  if (x[5] > 1.0) x[5] = 1.0;
153 
154  xOut = 1.;
155  simu->Fill();
156  if (i >= fMaxNofTrainSamples) break;
157  }
158  for (size_t i = 0; i < IM_wrong.size(); i++) {
159  x[0] = IM_wrong[i] / 0.1;
160  x[1] = OA_wrong[i] / 30;
161  x[2] = Angle_wrong[i] / 30;
162  x[3] = Z_wrong[i] / 100;
163  x[4] = Mom1_wrong[i] / 5;
164  x[5] = Mom2_wrong[i] / 5;
165 
166  if (x[0] > 1.0) x[0] = 1.0;
167  if (x[1] > 1.0) x[1] = 1.0;
168  if (x[2] > 1.0) x[2] = 1.0;
169  if (x[3] > 1.0) x[3] = 1.0;
170  if (x[4] > 1.0) x[4] = 1.0;
171  if (x[5] > 1.0) x[5] = 1.0;
172 
173  xOut = -1.;
174  simu->Fill();
175  if (i >= fMaxNofTrainSamples) break;
176  }
177 
178  TMultiLayerPerceptron network("x0,x1,x2,x3,x4,x5:12:xOut", simu, "Entry$+1");
179  network.Train(300, "text,update=10");
180  network.DumpWeights(
181  "../../../analysis/conversion2/KresAnalysis_ann_weights.txt");
182 
183 
184  Double_t params[6];
185  fNofWrongLikeCorrect = 0;
186  fNofCorrectLikeWrong = 0;
187 
188  for (size_t i = 0; i < IM_correct.size(); i++) {
189  params[0] = IM_correct[i] / 0.1;
190  params[1] = OA_correct[i] / 30;
191  params[2] = Angle_correct[i] / 30;
192  params[3] = Z_correct[i] / 100;
193  params[4] = Mom1_correct[i] / 5;
194  params[5] = Mom2_correct[i] / 5;
195 
196  if (params[0] > 1.0) params[0] = 1.0;
197  if (params[1] > 1.0) params[1] = 1.0;
198  if (params[2] > 1.0) params[2] = 1.0;
199  if (params[3] > 1.0) params[3] = 1.0;
200  if (params[4] > 1.0) params[4] = 1.0;
201  if (params[5] > 1.0) params[5] = 1.0;
202 
203  Double_t netEval = network.Evaluate(0, params);
204  fhAnnOutput_correct->Fill(netEval);
205  if (netEval < fAnnCut) fNofCorrectLikeWrong++;
206  }
207  for (size_t i = 0; i < IM_wrong.size(); i++) {
208  params[0] = IM_wrong[i] / 0.1;
209  params[1] = OA_wrong[i] / 30;
210  params[2] = Angle_wrong[i] / 30;
211  params[3] = Z_wrong[i] / 100;
212  params[4] = Mom1_wrong[i] / 5;
213  params[5] = Mom2_wrong[i] / 5;
214 
215  if (params[0] > 1.0) params[0] = 1.0;
216  if (params[1] > 1.0) params[1] = 1.0;
217  if (params[2] > 1.0) params[2] = 1.0;
218  if (params[3] > 1.0) params[3] = 1.0;
219  if (params[4] > 1.0) params[4] = 1.0;
220  if (params[5] > 1.0) params[5] = 1.0;
221 
222  Double_t netEval = network.Evaluate(0, params);
223  fhAnnOutput_wrong->Fill(netEval);
224  if (netEval >= fAnnCut) fNofWrongLikeCorrect++;
225  }
226 }
227 
228 
229 void CbmKresTrainAnn::Draw() {
230  cout << "nof correct pairs = " << IM_correct.size() << endl;
231  cout << "nof wrong pairs = " << IM_wrong.size() << endl;
232  cout << "wrong like correct = " << fNofWrongLikeCorrect
233  << ", wrong supp = " << (Double_t) IM_wrong.size() / fNofWrongLikeCorrect
234  << endl;
235  cout << "Correct like wrong = " << fNofCorrectLikeWrong
236  << ", correct lost eff = "
237  << 100. * (Double_t) fNofCorrectLikeWrong / IM_correct.size() << endl;
238 
239  Double_t cumProbFake = 0.;
240  Double_t cumProbTrue = 0.;
241  Int_t nofTrue = (Int_t) fhAnnOutput_correct->GetEntries();
242  Int_t nofFake = (Int_t) fhAnnOutput_wrong->GetEntries();
243 
244  for (Int_t i = 1; i <= fhAnnOutput_wrong->GetNbinsX(); i++) {
245  cumProbTrue += fhAnnOutput_correct->GetBinContent(i);
246  fhCumProb_correct->SetBinContent(i, 1. - (Double_t) cumProbTrue / nofTrue);
247 
248  cumProbFake += fhAnnOutput_wrong->GetBinContent(i);
249  fhCumProb_wrong->SetBinContent(i, (Double_t) cumProbFake / nofFake);
250  }
251 
252 
253  TCanvas* c1 =
254  new TCanvas("ann_correct_ann_output", "ann_correct_ann_output", 400, 400);
255  c1->SetTitle("ann_correct_ann_output");
256  fhAnnOutput_correct->Draw();
257 
258  TCanvas* c2 =
259  new TCanvas("ann_wrong_ann_output", "ann_wrong_ann_output", 400, 400);
260  c2->SetTitle("ann_wrong_ann_output");
261  fhAnnOutput_wrong->Draw();
262 
263  TCanvas* c3 =
264  new TCanvas("ann_correct_cum_prob", "ann_correct_cum_prob", 400, 400);
265  c3->SetTitle("ann_correct_cum_prob");
266  fhCumProb_correct->Draw();
267 
268  TCanvas* c4 =
269  new TCanvas("ann_wrong_cum_prob", "ann_wrong_cum_prob", 400, 400);
270  c4->SetTitle("ann_wrong_cum_prob");
271  fhCumProb_wrong->Draw();
272 }
273 
274 
275 void CbmKresTrainAnn::InitHistograms() {
276 
277  fhAnnOutput_correct = new TH1D(
278  "fhAnnOutput_correct", "ANN output;ANN output;Counter", 100, -1.2, 1.2);
279  fHists.push_back(fhAnnOutput_correct);
280  fhAnnOutput_wrong = new TH1D(
281  "fhAnnOutput_wrong", "ANN output;ANN output;Counter", 100, -1.2, 1.2);
282  fHists.push_back(fhAnnOutput_wrong);
283 
284  fhCumProb_correct = new TH1D("fhCumProb_correct",
285  "ANN output;ANN output;Cumulative probability",
286  100,
287  -1.2,
288  1.2);
289  fHists.push_back(fhCumProb_correct);
290  fhCumProb_wrong = new TH1D("fhCumProb_wrong",
291  "ANN output;ANN output;Cumulative probability",
292  100,
293  -1.2,
294  1.2);
295  fHists.push_back(fhCumProb_wrong);
296 }
CbmKresTrainAnn::fMaxNofTrainSamples
unsigned int fMaxNofTrainSamples
Definition: CbmKresTrainAnn.h:37
CbmKresTrainAnn::Mom2_wrong
vector< double > Mom2_wrong
Definition: CbmKresTrainAnn.h:53
CbmKresTrainAnn::TrainAndTestAnn
void TrainAndTestAnn()
Definition: CbmKresTrainAnn.cxx:125
CbmKresTrainAnn::fAnnCut
double fAnnCut
Definition: CbmKresTrainAnn.h:38
CbmKresTrainAnn::OA_correct
vector< double > OA_correct
Definition: CbmKresTrainAnn.h:43
CbmKresTrainAnn::fHists
vector< TH1 * > fHists
Definition: CbmKresTrainAnn.h:55
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmKresTrainAnn::Mom1_correct
vector< double > Mom1_correct
Definition: CbmKresTrainAnn.h:46
CbmKresTrainAnn::Angle_wrong
vector< double > Angle_wrong
Definition: CbmKresTrainAnn.h:50
CbmKresTrainAnn::IM_wrong
vector< double > IM_wrong
Definition: CbmKresTrainAnn.h:48
CbmKresTrainAnn::Z_correct
vector< double > Z_correct
Definition: CbmKresTrainAnn.h:45
CbmKresTrainAnn::IM_correct
vector< double > IM_correct
Definition: CbmKresTrainAnn.h:42
CbmDrawHist.h
Helper functions for drawing 1D and 2D histograms and graphs.
CbmKresTrainAnn::CbmKresTrainAnn
CbmKresTrainAnn()
Definition: CbmKresTrainAnn.cxx:40
CbmKresTrainAnn::Angle_correct
vector< double > Angle_correct
Definition: CbmKresTrainAnn.h:44
CbmKresTrainAnn::InitHistograms
void InitHistograms()
Definition: CbmKresTrainAnn.cxx:275
CbmKresTrainAnn::Mom1_wrong
vector< double > Mom1_wrong
Definition: CbmKresTrainAnn.h:52
CbmKresTrainAnn::fhCumProb_wrong
TH1D * fhCumProb_wrong
Definition: CbmKresTrainAnn.h:59
CbmKresTrainAnn::fhAnnOutput_wrong
TH1D * fhAnnOutput_wrong
Definition: CbmKresTrainAnn.h:57
CbmKresTrainAnn.h
CbmKresTrainAnn::fNofWrongLikeCorrect
int fNofWrongLikeCorrect
Definition: CbmKresTrainAnn.h:39
CbmKresTrainAnn::~CbmKresTrainAnn
virtual ~CbmKresTrainAnn()
Definition: CbmKresTrainAnn.cxx:63
CbmKresTrainAnn::Exec
void Exec(int event, int IdForANN, double InvariantMass, double OpeningAngle, double PlaneAngle_last, double ZPos, TVector3 Momentum1, TVector3 Momentum2)
Definition: CbmKresTrainAnn.cxx:67
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmKresTrainAnn::Z_wrong
vector< double > Z_wrong
Definition: CbmKresTrainAnn.h:51
CbmKresTrainAnn::fhCumProb_correct
TH1D * fhCumProb_correct
Definition: CbmKresTrainAnn.h:58
CbmKresTrainAnn::OA_wrong
vector< double > OA_wrong
Definition: CbmKresTrainAnn.h:49
CbmKresTrainAnn::fhAnnOutput_correct
TH1D * fhAnnOutput_correct
Definition: CbmKresTrainAnn.h:56
CbmKresTrainAnn::Init
void Init()
Definition: CbmKresTrainAnn.cxx:65
CbmKresTrainAnn::fNofCorrectLikeWrong
int fNofCorrectLikeWrong
Definition: CbmKresTrainAnn.h:40
CbmKresTrainAnn::Draw
void Draw()
Definition: CbmKresTrainAnn.cxx:229
CbmKresTrainAnn::Mom2_correct
vector< double > Mom2_correct
Definition: CbmKresTrainAnn.h:47