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