ATLAS Offline Software
positions/trainNN.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include <TTree.h>
6 #include <TFile.h>
7 #include <TCanvas.h>
8 #include "TChain.h"
9 #include <TH1F.h>
10 #include <TLegend.h>
11 #include <TPad.h>
12 
13 #include "../TJetNet.h"
14 #include "../doNormalization.C"
15 #include "Riostream.h"
16 #include "../TNetworkToHistoTool.h"
17 #include "TSystem.h"
18 
19 #include "../TTrainedNetwork.h"
20 
21 #include "TMatrixD.h"
22 #include "TVectorD.h"
23 #include "trainNN.h"
24 
25 #include <vector>
26 #include <utility>
27 #include <algorithm>
28 #include <string>
29 #include <cmath>
30 #include <iostream>
31 #include <stdexcept>
32 
33 using namespace std;
34 
35 Double_t sigmoid(Double_t x)
36 {
37  return 1./(1.+exp(-2*x));
38 }
39 
40 using namespace std;
41 bool isBadCluster( int sizeX, int nParticles ){
42 
43 
44  if(sizeX==-100){return true;}
45  if(nParticles==0){return true;}
46  if(nParticles!=1 && nParticles!=2 && nParticles!=3){return true;}
47 
48 
49 
50  return false;
51 
52 
53 }
54 
55 
56 bool skipSingle(int nParticles,int iClus, int dilutionFactor){
57 
58 
59  // this reduces the number of single particle clusters of a factor k
60  int k = 60;
61 
62 
63  if (nParticles==1)
64  {
65  if ((iClus %(dilutionFactor*k) != 0) && (iClus % (dilutionFactor*k) != 1))
66  {
67  return true;
68  }
69  }
70 
71  return false;
72 
73 }
74 
75 
76 
77 bool badTrackInfo(bool useTrackEstimate, double theta){
78 
79  if(useTrackEstimate && theta==0){return true;}
80  return false;
81 
82 }
83 
84 
85 
86 int main()
87 {
88  trainNN(
89  "../reduceddatasets/reduceddataset_Cone4H1TopoParticleJets_forNN.root",
90  "dummy",
91  10000,
92  200,
93  10,
94  10,
95  0);
96  return 0;
97 }
98 
99 
100 
101 void trainNN(TString inputfile,
102  TString outputclass,
103  int nIterations,
104  int dilutionFactor,
105  int nodesFirstLayer,
106  int nodesSecondLayer,
107  int restartTrainingFrom,
108  int nParticlesTraining,
109  bool useTrackEstimate,
110  int nPatternsPerUpdate,
111  double learningRate,
112  double learningRateDecrease,
113  double learningRateMomentum) {
114 
115  double bweight=1;
116  double cweight=1.;
117  double lweight=1;
118 
119  gROOT->SetStyle("Plain");
120 
121  cout << "starting with settings: " << endl;
122  cout << " nIterations: " << nIterations << endl;
123  cout << " dilutionFactor: " << dilutionFactor << endl;
124  cout << " nodesFirstLayer: " << nodesFirstLayer << endl;
125  cout << " nodesSecondLayer: " << nodesSecondLayer << endl;
126 
127 // TFile *file = TFile::Open(inputfile);
128 // TTree *simu = (TTree*)file->Get("Validation/NNinput");
129 
130  TChain *myChain = new TChain("Validation/NNinput");
131 
132 
133 if(!useTrackEstimate){
134  #include "../files.txt"
135 }
136 
137 if(useTrackEstimate){
138  #include "../filesOnTrack.txt"
139 }
140  TChain* simu=myChain;
141 
142  std::cout << " Training sample obtained... " << std::endl;
143 
144  vector<int> *NN_sizeX;
145  vector<int> *NN_sizeY;
146  vector<vector<float> > *NN_matrixOfToT;
147  vector<vector<float> > *NN_vectorOfPitchesY;
148  vector<int> *NN_ClusterPixLayer;
149  vector<int> *NN_ClusterPixBarrelEC;
150  vector<float> *NN_phiBS;
151  vector<float> *NN_thetaBS;
152  vector<float> *NN_etaModule;
153  vector<bool> *NN_useTrackInfo;
154  vector<int> *NN_columnWeightedPosition;
155  vector<int> *NN_rowWeightedPosition;
156  vector<vector<float> > *NN_positionX;
157  vector<vector<float> > *NN_positionY;
158  vector<vector<float> > *NN_position_idX;
159  vector<vector<float> > *NN_position_idY;
160  vector<vector<float> > *NN_theta;
161  vector<vector<float> > *NN_phi;
162 
163  // List of branches
164  TBranch *b_NN_sizeX;
165  TBranch *b_NN_sizeY;
166  TBranch *b_NN_matrixOfToT;
167  TBranch *b_NN_vectorOfPitchesY;
168  TBranch *b_NN_ClusterPixLayer;
169  TBranch *b_NN_ClusterPixBarrelEC;
170  TBranch *b_NN_phiBS;
171  TBranch *b_NN_thetaBS;
172  TBranch *b_NN_etaModule;
173  TBranch *b_NN_useTrackInfo;
174  TBranch *b_NN_columnWeightedPosition;
175  TBranch *b_NN_rowWeightedPosition;
176  TBranch *b_NN_positionX;
177  TBranch *b_NN_positionY;
178  TBranch *b_NN_position_idX;
179  TBranch *b_NN_position_idY;
180  TBranch *b_NN_theta;
181  TBranch *b_NN_phi;
182 
183 
184 
185  NN_sizeX = 0;
186  NN_sizeY = 0;
187  NN_matrixOfToT = 0;
188  NN_vectorOfPitchesY = 0;
189  NN_ClusterPixLayer = 0;
190  NN_ClusterPixBarrelEC = 0;
191  NN_phiBS = 0;
192  NN_thetaBS = 0;
193  NN_etaModule = 0;
194  NN_useTrackInfo = 0;
195  NN_columnWeightedPosition = 0;
196  NN_rowWeightedPosition = 0;
197  NN_positionX = 0;
198  NN_positionY = 0;
199  NN_position_idX = 0;
200  NN_position_idY = 0;
201  NN_theta = 0;
202  NN_phi = 0;
203  // Set branch addresses and branch pointers
204  // if (!tree) return 0;
205  // TTree* simu = tree;
206  // fCurrent = -1;
207  simu->SetMakeClass(1);
208 
209  simu->SetBranchAddress("NN_sizeX", &NN_sizeX, &b_NN_sizeX);
210  simu->SetBranchAddress("NN_sizeY", &NN_sizeY, &b_NN_sizeY);
211  simu->SetBranchAddress("NN_matrixOfToT", &NN_matrixOfToT, &b_NN_matrixOfToT);
212  simu->SetBranchAddress("NN_vectorOfPitchesY", &NN_vectorOfPitchesY, &b_NN_vectorOfPitchesY);
213  simu->SetBranchAddress("NN_ClusterPixLayer", &NN_ClusterPixLayer, &b_NN_ClusterPixLayer);
214  simu->SetBranchAddress("NN_ClusterPixBarrelEC", &NN_ClusterPixBarrelEC, &b_NN_ClusterPixBarrelEC);
215  simu->SetBranchAddress("NN_phiBS", &NN_phiBS, &b_NN_phiBS);
216  simu->SetBranchAddress("NN_thetaBS", &NN_thetaBS, &b_NN_thetaBS);
217  simu->SetBranchAddress("NN_etaModule", &NN_etaModule, &b_NN_etaModule);
218  simu->SetBranchAddress("NN_useTrackInfo", &NN_useTrackInfo, &b_NN_useTrackInfo);
219  simu->SetBranchAddress("NN_columnWeightedPosition", &NN_columnWeightedPosition, &b_NN_columnWeightedPosition);
220  simu->SetBranchAddress("NN_rowWeightedPosition", &NN_rowWeightedPosition, &b_NN_rowWeightedPosition);
221  simu->SetBranchAddress("NN_positionX", &NN_positionX, &b_NN_positionX);
222  simu->SetBranchAddress("NN_positionY", &NN_positionY, &b_NN_positionY);
223  simu->SetBranchAddress("NN_position_idX", &NN_position_idX, &b_NN_position_idX);
224  simu->SetBranchAddress("NN_position_idY", &NN_position_idY, &b_NN_position_idY);
225  simu->SetBranchAddress("NN_theta", &NN_theta, &b_NN_theta);
226  simu->SetBranchAddress("NN_phi", &NN_phi, &b_NN_phi);
227 
228 
229  cout << "Branches set..." << endl;
230 
231  TString filterTrain("Entry$%");
232  filterTrain+=dilutionFactor;
233  filterTrain+="==0";
234 
235  TString filterTest("Entry$%");
236  filterTest+=dilutionFactor;
237  filterTest+="==1";
238 
239  int* nneurons;
240  int nlayer=3;
241 
242  simu->GetEntry(0);
243 
244 
245  cout << "First entry..." << endl;
246 
247  Int_t sizeX=-7;
248  Int_t sizeY=-7;
249 
250 
251  // loop over the clusters loking for the first cluster properly saved
252  for( unsigned int clus =0; clus<NN_sizeX->size(); clus++ ){
253 
254  sizeX = (*NN_sizeX)[clus];
255  sizeY = (*NN_sizeY)[clus];
256 
257  if(sizeX>0)break;
258 
259  }
260 
261  cout << "Size obtained" << endl;
262 
263 
264  int numberinputs=sizeX*(sizeY+1)+4;
265  if (!useTrackEstimate)
266  {
267  numberinputs=sizeX*(sizeY+1)+5;
268  }
269 
270  int numberoutputs=2*nParticlesTraining;
271 
272  if (nodesSecondLayer!=0)
273  {
274  nlayer=4;
275  }
276 
277  if (nodesSecondLayer!=0)
278  {
279  nneurons=new int[4];
280  }
281  else
282  {
283  nneurons=new int[3];
284  }
285 
286  nneurons[0]=numberinputs;
287 
288  nneurons[1]=nodesFirstLayer;
289 
290  if (nodesSecondLayer!=0)
291  {
292  nneurons[2]=nodesSecondLayer;
293  nneurons[3]=numberoutputs;//number of output nodes
294  }
295  else
296  {
297  nneurons[2]=numberoutputs;//number of output nodes
298  }
299 
300  for (int i=0;i<nlayer;i++)
301  {
302  cout << " layer i: " << i << " number neurons: " << nneurons[i] << endl;
303  }
304 
305 
306  // float eventWeight(0);
307  float trainingError(0);
308  float testError(0);
309 
310  //setting learning parameters
311 
312  cout << " now providing training events " << endl;
313 
314  Int_t numberTrainingEvents=0;
315  Int_t numberTestingEvents=0;
316 
317  int iClus=0;
318  int part_0=0;
319  int part_1=0;
320  int part_2=0;
321  int part_3=0;
322 
323  int totalN=simu->GetEntries();
324 
325 // totalN=10000;
326 
327  // Loop over entries:
328  for (Int_t i = 0; i < totalN; i++) {
329 
330  if (i % 50000 == 0 ) {
331  std::cout << " Counting training / testing events in sample. Looping over event " << i << std::endl;
332  }
333 
334  simu->GetEntry(i);
335 
336  for( unsigned int clus =0; clus<NN_sizeX->size(); clus++ ){
337 
338  vector<float> *matrixOfToT=0;
339  vector<float> *vectorOfPitchesY=0;
340 
341  Float_t phiBS;
342  Float_t thetaBS;
343  Float_t etaModule;
344  Int_t ClusterPixLayer;
345  Int_t ClusterPixBarrelEC;
346 
347  std::vector<float> * positionX=0;
348  std::vector<float> * positionY=0;
349  std::vector<float> * thetaTr=0;
350  std::vector<float> * phiTr=0;
351 
352  std::vector<float> positionX_reorder;
353  std::vector<float> positionY_reorder;
354 
355 
356  sizeX = (*NN_sizeX)[clus];
357  positionX =&(*NN_positionX)[clus];
358  int nParticles = positionX->size();
359 
360  thetaTr = &(*NN_theta)[clus];
361  phiTr = &(*NN_phi)[clus];
362 
363  if (nParticlesTraining!=nParticles)
364  {
365  continue;
366  }
367  if (isBadCluster(sizeX, nParticles ) )continue;
368 
369  // loop over the particles;
370  for( unsigned int P = 0; P < positionX->size(); P++){
371  double theta = (*thetaTr)[P];
372  if (theta!=theta) continue;
373 
374  iClus++;
375 
376  if ( badTrackInfo(useTrackEstimate, theta ) )continue;
377 
378 
379 
380  if (iClus%dilutionFactor==0) numberTrainingEvents+=1;
381  if (iClus%dilutionFactor==1) numberTestingEvents+=1;
382 
383  if (iClus%dilutionFactor==1 && nParticles==1 ) part_1++;
384  if (iClus%dilutionFactor==1 && nParticles==2 ) part_2++;
385  if (iClus%dilutionFactor==1 && nParticles==3 ) part_3++;
386 
387 
388 
389  }// end loop over th particles
390  }// end loop over cluster
391  }// end Loop over entries
392 
393 
394 
395  cout << " N. training events: " << numberTrainingEvents <<
396  " N. testing events: " << numberTestingEvents << endl;
397 
398  // if(numberTrainingEvents!=numberTestingEvents)return;
399 
400  cout << "now start to setup the network..." << endl;
401 
402  TJetNet* jn = new TJetNet( numberTestingEvents, numberTrainingEvents, nlayer, nneurons );
403  // return;
404  cout << " setting learning method... " << endl;
405 
406  // jn->SetMSTJN(4,12); Fletscher-Rieves (Scaled Conj Grad)
407 
408 
409 // int nPatternsPerUpdate=1;
410 
411  jn->SetPatternsPerUpdate( nPatternsPerUpdate );
412  jn->SetUpdatesPerEpoch( (int)std::floor((float)numberTrainingEvents/(float)nPatternsPerUpdate) );
413  jn->SetUpdatingProcedure( 0 );
414  jn->SetErrorMeasure( 0 );
415  jn->SetActivationFunction( 1 );
416 
417  jn->SetLearningRate( learningRate );//0.8
418 // jn->SetLearningRate( 1 );//0.8
419 // jn->SetMomentum( 0.2 );//0.3 //is now 0.5
420  jn->SetMomentum( learningRateMomentum );//0.3 //is now 0.5
421  jn->SetInitialWeightsWidth( 1. );
422  // jn->SetLearningRateDecrease( 0.992 );
423 
424  jn->SetLearningRateDecrease( learningRateDecrease );//0.992
425 // jn->SetLearningRateDecrease( 0.995 );//0.992
426  jn->SetActivationFunction(4,nlayer-1-1);
427 
428 
429 
430  cout << " setting pattern for training events " << endl;
431 
432  int trainSampleNumber=0;
433  int testSampleNumber=1;
434 
435  cout << " copying over training events " << endl;
436 
437  int counter=0;
438  int counter0=0;
439  int counter1=0;
440 
441  //input and output of first event
442  vector<double> inputVar;
443  vector<double> outputVar;
444 
445 
446  iClus=0;
447  for (Int_t i = 0; i < totalN; i++) {
448 
449  if (i % 100000 == 0 ) {
450  std::cout << " Copying over training events. Looping over event " << i << std::endl;
451  }
452 
453 
454  simu->GetEntry(i);
455 
456  for( unsigned int clus =0; clus<NN_sizeX->size(); clus++ ){
457 
458  vector<float> *matrixOfToT=0;
459  vector<float> *vectorOfPitchesY=0;
460 
461  Float_t phiBS;
462  Float_t thetaBS;
463  Float_t etaModule;
464  Int_t ClusterPixLayer;
465  Int_t ClusterPixBarrelEC;
466 
467  std::vector<float> * position_idX=0;
468  std::vector<float> * position_idY=0;
469  std::vector<float> * thetaTr=0;
470  std::vector<float> * phiTr=0;
471 
472 
473  sizeX = (*NN_sizeX)[clus];
474 
475  sizeY = (*NN_sizeY)[clus];
476  matrixOfToT=&(*NN_matrixOfToT)[clus];
477  vectorOfPitchesY=&(*NN_vectorOfPitchesY)[clus];
478 
479  phiBS = (*NN_phiBS)[clus];
480  thetaBS =(*NN_thetaBS)[clus];
481  etaModule =(*NN_etaModule)[clus];
482 
483  ClusterPixLayer=(*NN_ClusterPixLayer)[clus];
484  ClusterPixBarrelEC = (*NN_ClusterPixBarrelEC)[clus];
485 
486  position_idX =&(*NN_position_idX)[clus];
487  position_idY =&(*NN_position_idY)[clus];
488 
489  int nParticles = position_idX->size();
490 
491  if (nParticlesTraining!=nParticles)
492  {
493  continue;
494  }
495 
496  thetaTr = &(*NN_theta)[clus];
497  phiTr = &(*NN_phi)[clus];
498 
499 
500  if(isBadCluster(sizeX, nParticles ) )continue;
501 
502  for( unsigned int P = 0; P < position_idX->size(); P++){
503 
504  double theta = (*thetaTr)[P];
505  double phi = (*phiTr)[P];
506  if (theta!=theta) continue;
507 
508  if (ClusterPixBarrelEC==2)
509  {
510  theta=-theta;
511  phi=-phi;
512  }
513 
514  iClus++;
515 
516 
517 
518  if ( badTrackInfo(useTrackEstimate, theta ) )continue;
519 
520 
521 
522  if (matrixOfToT->size()!=sizeX*sizeY)
523  {
524  std::cout << " Event: " << i << " PROBLEM: size Y is: " << matrixOfToT->size() << std::endl;
525  throw std::runtime_error("Error in positions/trainNN.cxx");
526  }
527 
528  // loop over elements of matrixOfTot which is actually a vector
529  for(unsigned int ME =0; ME < matrixOfToT->size(); ME++){
530 
531  if (iClus%dilutionFactor==0) jn->SetInputTrainSet( counter0, ME, norm_ToT((*matrixOfToT)[ME]));
532  if (iClus%dilutionFactor==1) jn->SetInputTestSet( counter1, ME, norm_ToT((*matrixOfToT)[ME]));
533 
534  if ((*matrixOfToT)[ME] != (*matrixOfToT)[ME])
535  {
536  cout << "ME n. " << ME << " is: " << (*matrixOfToT)[ME] << endl;
537  throw std::runtime_error("Error in positions/trainNN.cxx");
538  }
539 
540  if (counter0 == 0) std::cout << " element: " << ME << " ToT set to: " << norm_ToT((*matrixOfToT)[ME]) << std::endl;
541  if (iClus%dilutionFactor==1&&counter1==0) inputVar.push_back(norm_ToT((*matrixOfToT)[ME]));
542 
543  }
544 
545 
546  // loop over vector of pitches
547  for (int s=0;s<sizeY;s++)
548  {
549  if (iClus%dilutionFactor==0) jn->SetInputTrainSet( counter0, sizeX*sizeY+s, norm_pitch((*vectorOfPitchesY)[s]));
550  if (iClus%dilutionFactor==1) jn->SetInputTestSet( counter1, sizeX*sizeY+s, norm_pitch((*vectorOfPitchesY)[s]));
551 
552  if ((*vectorOfPitchesY)[s]!=(*vectorOfPitchesY)[s])
553  {
554  cout << " pitchY: " << (*vectorOfPitchesY)[s] << endl;
555  throw std::runtime_error("Error in positions/trainNN.cxx");;
556  }
557  if (counter0 == 0) std::cout << " s: " << s << " pitch set to: " << norm_pitch((*vectorOfPitchesY)[s]) << std::endl;
558 
559  if (iClus%dilutionFactor==1&&counter1==0) inputVar.push_back(norm_pitch((*vectorOfPitchesY)[s]));
560  }
561 
562  if (iClus%dilutionFactor==0) jn->SetInputTrainSet( counter0, (sizeX+1)*sizeY, norm_layerNumber(ClusterPixLayer));
563  if (iClus%dilutionFactor==0) jn->SetInputTrainSet( counter0, (sizeX+1)*sizeY+1, norm_layerType(ClusterPixBarrelEC));
564 
565  if (ClusterPixLayer!=ClusterPixLayer || ClusterPixBarrelEC!=ClusterPixBarrelEC)
566  {
567  cout << " ClusterPixLayer: " << ClusterPixLayer << " ClusterPixBarrelEC " << ClusterPixBarrelEC << endl;
568  throw std::runtime_error("Error in positions/trainNN.cxx");
569  }
570 
571  if (iClus%dilutionFactor==1) jn->SetInputTestSet( counter1, (sizeX+1)*sizeY, norm_layerNumber(ClusterPixLayer));
572  if (iClus%dilutionFactor==1) jn->SetInputTestSet( counter1, (sizeX+1)*sizeY+1, norm_layerType(ClusterPixBarrelEC));
573 
574  if (iClus%dilutionFactor==1&&counter1==0) {
575  inputVar.push_back(norm_layerNumber(ClusterPixLayer));
576  inputVar.push_back(norm_layerType(ClusterPixBarrelEC));
577  }
578 
579 
580  if (counter0 == 0) std::cout << " ClusterPixLayer " << norm_layerNumber(ClusterPixLayer) << " ClusterPixBarrelEC " << norm_layerType(ClusterPixBarrelEC) << std::endl;
581 
582  if (useTrackEstimate)
583  {
584  if (iClus%dilutionFactor==0) jn->SetInputTrainSet( counter0, (sizeX+1)*sizeY+2, norm_phi(phi) );
585  if (iClus%dilutionFactor==0) jn->SetInputTrainSet( counter0, (sizeX+1)*sizeY+3, norm_theta(theta) );
586 
587  if (phi!=phi) {
588  cout << " phi: " << phi << endl;
589  throw std::runtime_error("Error in positions/trainNN.cxx");
590  }
591 
592  if (theta!=theta) {
593  cout << " theta: " << theta << endl;
594  throw std::runtime_error("Error in positions/trainNN.cxx");
595  }
596 
597 
598 
599  if (iClus%dilutionFactor==1) jn->SetInputTestSet( counter1, (sizeX+1)*sizeY+2, norm_phi(phi) );
600  if (iClus%dilutionFactor==1) jn->SetInputTestSet( counter1, (sizeX+1)*sizeY+3, norm_theta(theta) );
601 
602  if (counter0==0) std::cout << " phi " << norm_phi(phi) << " theta: " << norm_theta(theta) << std::endl;
603 
604  if (iClus%dilutionFactor==1&&counter1==0) {
605  inputVar.push_back(norm_phi(phi));
606  inputVar.push_back(norm_theta(theta));
607  }
608  }
609  else
610  {
611  if (iClus%dilutionFactor==0) jn->SetInputTrainSet( counter0, (sizeX+1)*sizeY+2, norm_phiBS(phiBS) );
612  if (iClus%dilutionFactor==0) jn->SetInputTrainSet( counter0, (sizeX+1)*sizeY+3, norm_thetaBS(thetaBS) );
613  if (iClus%dilutionFactor==0) jn->SetInputTrainSet( counter0, (sizeX+1)*sizeY+4, norm_etaModule(etaModule) );
614 
615  if (iClus%dilutionFactor==1) jn->SetInputTestSet( counter1, (sizeX+1)*sizeY+2, norm_phiBS(phiBS) );
616  if (iClus%dilutionFactor==1) jn->SetInputTestSet( counter1, (sizeX+1)*sizeY+3, norm_thetaBS(thetaBS) );
617  if (iClus%dilutionFactor==1) jn->SetInputTestSet( counter1, (sizeX+1)*sizeY+4, norm_etaModule(etaModule) );
618 
619 
620  if (iClus%dilutionFactor==1&&counter1==0) {
621  inputVar.push_back(norm_phiBS(phiBS));
622  inputVar.push_back(norm_thetaBS(thetaBS));
623  inputVar.push_back(norm_etaModule(etaModule) );
624  }
625 
626  if (counter0==0) std::cout <<
627  " phiBS " << norm_phiBS(phiBS) <<
628  " thetaBS: " << norm_thetaBS(thetaBS) <<
629  " etaModule: " << norm_etaModule(etaModule) << std::endl;
630  }
631 
632 
633  vector<float> xPositions = *position_idX;
634 
635  // for (int e=0;e<nParticles;e++)
636  // {
637  // xPositions.push_back((*positionsX)[e]);
638  // }
639 
640  std::sort(xPositions.begin(),xPositions.end());
641 
642  for (int o=0;o<xPositions.size();o++)
643  {
644 
645  if (iClus%dilutionFactor==0) jn->SetOutputTrainSet(counter0, 2*o, norm_posX(xPositions[o]));
646  if (iClus%dilutionFactor==1) jn->SetOutputTestSet(counter1, 2*o, norm_posX(xPositions[o]));
647 
648  if (xPositions[o]!=xPositions[o])
649  {
650  cout << "pos: " << xPositions[o] << endl;
651  throw std::runtime_error("Error in positions/trainNN.cxx");
652  }
653 
654  if (counter0==0) std::cout << " output node: " << 2*o << " set to: " << norm_posX(xPositions[o]) << endl;
655 
656  if (iClus%dilutionFactor==1&&counter1==0) { outputVar.push_back(norm_posX(xPositions[o]));}
657 
658  double corry=-1000;
659 
660  for (int e=0;e<nParticles;e++)
661  {
662  if (fabs((*position_idX)[e]-xPositions[o])<1e-10)
663  {
664  if (fabs(corry+1000)>1e-6)
665  {
666  cout << " Value find more than once! " << endl;
667  for (int p=0;p<xPositions.size();p++)
668  {
669  cout << " X n. : " << p << " is: " << xPositions[p] << endl;
670  }
671  }
672  corry=(*position_idY)[e];
673  }
674  }
675 
676  if (fabs(corry+1000)<1e-6) {
677  cout << " could not find original X pos. " << endl;
678  throw std::runtime_error("Error in positions/trainNN.cxx");
679  }
680 
681  if (iClus%dilutionFactor==0) jn->SetOutputTrainSet(counter0, 2*o+1, norm_posY(corry));
682  if (iClus%dilutionFactor==1) jn->SetOutputTestSet(counter1, 2*o+1, norm_posY(corry));
683 
684  if (corry!=corry) {
685  cout << " posY " << corry << endl;
686  throw std::runtime_error("Error in positions/trainNN.cxx");
687  }
688 
689 
690  if (counter0==0) std::cout << " output node: " << 2*o+1 << " set to: " << norm_posY(corry) << endl;
691  if (iClus%dilutionFactor==1&&counter1==0) { outputVar.push_back(norm_posY(corry));}
692 
693  }
694 
695  if (iClus%dilutionFactor==0) jn->SetEventWeightTrainSet( counter0, 1 );
696  if (iClus%dilutionFactor==1) jn->SetEventWeightTestSet( counter1, 1 );
697 
698 
699  if (iClus%dilutionFactor==0){counter0+=1;}
700  if (iClus%dilutionFactor==1){counter1+=1;}
701 
702 
703 
704  }
705  }
706  }
707 
708  if (counter0!=numberTrainingEvents)
709  {
710  cout << " counter up to: " << counter0 << " while events in training sample are " << numberTrainingEvents << endl;
711  return;
712  }
713 
714  if (counter1!=numberTestingEvents)
715  {
716  cout << " counter up to: " << counter1 << " while events in testing sample are " << numberTestingEvents << endl;
717  return;
718  }
719 
720  // return;
721 
722  jn->Shuffle(true,false);
723 
724  if (restartTrainingFrom==0)
725  {
726  jn->Init();
727  // jn->DumpToFile("WeightsInitial.w");
728  }
729  else
730  {
731  TString name("Weights");
732  name+=restartTrainingFrom;
733  name+=".w";
734 
735  jn->ReadFromFile(name);
736  }
737 
738 
739 
740  float minimumError=1e10;
741  int epochesWithRisingError=0;
742  int epochWithMinimum=0;
743 
744  int updatesPerEpoch=jn->GetUpdatesPerEpoch();
745 
746  //prepare output stream
747 
748  TString directory("weights");
749  directory+="_nPat";
750  directory+=nPatternsPerUpdate;
751  directory+="_rate";
752  directory+=(int)(learningRate*100);
753  directory+="_rdec";
754  directory+=(int)(100*learningRateDecrease);
755  directory+="_mom";
756  directory+=(int)(100*learningRateMomentum);
757  directory+="_";
758  directory+=nParticlesTraining;
759  if (useTrackEstimate)
760  {
761  directory+="_withTracks";
762  }
763 
764  TString command("mkdir ");
765  command+=directory;
766 
767  gSystem->Exec(command);
768 
769  TString nameCronology=directory;
770 
771  nameCronology+="/trainingCronology.txt";
772 
773  ofstream cronology(nameCronology,ios_base::out);//|ios_base::app);
774 
775  cronology << "-------------SETTINGS----------------" << endl;
776  cronology << "Epochs: " << jn->GetEpochs() << std::endl;
777  cronology << "Updates Per Epoch: " << jn->GetUpdatesPerEpoch() << std::endl;
778  cronology << "Updating Procedure: " << jn->GetUpdatingProcedure() << std::endl;
779  cronology << "Error Measure: " << jn->GetErrorMeasure() << std::endl;
780  cronology << "Patterns Per Update: " << jn->GetPatternsPerUpdate() << std::endl;
781  cronology << "Learning Rate: " << jn->GetLearningRate() << std::endl;
782  cronology << "Momentum: " << jn->GetMomentum() << std::endl;
783  cronology << "Initial Weights Width: " << jn->GetInitialWeightsWidth() << std::endl;
784  cronology << "Learning Rate Decrease: " << jn->GetLearningRateDecrease() << std::endl;
785  cronology << "Activation Function: " << jn->GetActivationFunction() << std::endl;
786  cronology << "-------------LAYOUT------------------" << endl;
787  cronology << "Input variables: " << jn->GetInputDim() << endl;
788  cronology << "Output variables: " << jn->GetOutputDim() << endl;
789  cronology << "Hidden layers: " << jn->GetHiddenLayerDim() << endl;
790  cronology << "Layout : ";
791  for (Int_t s=0;s<jn->GetHiddenLayerDim()+2;++s)
792  {
793  cronology << jn->GetHiddenLayerSize(s);
794  if (s<jn->GetHiddenLayerDim()+1) cronology << "-";
795  }
796  cronology << endl;
797  cronology << "--------------HISTORY-----------------" << endl;
798  cronology << "History of iterations: " << endl;
799  cronology.close();
800 
801  //prepare training histo
802  TH1F* histoTraining=new TH1F("training","training",(int)std::floor((float)nIterations/10.+0.5),1,std::floor((float)nIterations/10.+1.5));
803  TH1F* histoTesting=new TH1F("testing","testing",(int)std::floor((float)nIterations/10.+0.5),1,std::floor((float)nIterations/10.+1.5));
804 
805  double maximumTrain=0;
806  double minimumTrain=1e10;
807 
808  for(int epoch=restartTrainingFrom+1;epoch<=nIterations;++epoch)
809  {
810  if (epoch!=restartTrainingFrom+1)
811  {
812  trainingError = jn->Train();
813  }
814 
815  if (epoch%10==0 || epoch==restartTrainingFrom+1)
816  {
817 
818  cronology.open(nameCronology,ios_base::app);
819 
820  testError = jn->Test();
821 
822  if (trainingError>maximumTrain) maximumTrain=trainingError;
823  if (testError>maximumTrain) maximumTrain=testError;
824  if (trainingError<minimumTrain) minimumTrain=trainingError;
825  if (testError<minimumTrain) minimumTrain=testError;
826 
827 
828  histoTraining->Fill(epoch/10.,trainingError);
829  histoTesting->Fill(epoch/10.,testError);
830 
831  if (testError<minimumError)
832  {
833  minimumError=testError;
834  epochesWithRisingError=0;
835  epochWithMinimum=epoch;
836  }
837  else
838  {
839  epochesWithRisingError+=10;
840 
841  }
842 
843 
844  if (epochesWithRisingError>300)
845  {
846  if (trainingError<minimumError)
847  {
848  cout << " End of training. Minimum already on epoch: " << epochWithMinimum << endl;
849  cronology << " End of training. Minimum already on epoch: " << epochWithMinimum << endl;
850  break;
851  }
852  }
853 
854  cronology << "Epoch: [" << epoch <<
855  "] Error: " << trainingError <<
856  " Test: " << testError << endl;
857 
858  cout << "Epoch: [" << epoch <<
859  "] Error: " << trainingError <<
860  " Test: " << testError << endl;
861 
862  cronology.close();
863 
864  TString name=directory;
865  name+="/Weights";
866  name+=epoch;
867  name+=".root";
868 
869  TFile* file=new TFile(name,"recreate");
870  TTrainedNetwork* trainedNetwork=jn->createTrainedNetwork();
871 
872 
873  jn->Evaluate( 0 ); //evaluate the first test pattern
874  for (int z=0;z<nParticlesTraining;z++)
875  {
876  std::cout << "output "<<z<<" x: " << jn->GetOutput(2*z) << " y: " << jn->GetOutput(2*z+1) << endl;
877  }
878 
879  std::vector<double> myTestOutput=trainedNetwork->calculateOutputValues(inputVar);
880  for (int z=0;z<nParticlesTraining;z++)
881  {
882  std::cout << "output TTNet "<<z<<" x: " << myTestOutput[2*z] <<
883  " y: " << myTestOutput[2*z+1] << endl;
884  }
885  for (int z=0;z<nParticlesTraining;z++)
886  {
887  std::cout << "should be " << z << "x: " << outputVar[2*z] << " y: " << outputVar[2*z+1] << endl;
888  }
889 
890  trainedNetwork->Write();
891  file->Write();
892  file->Close();
893  delete file;
894 
895  /*
896  TFile* file2=new TFile(name);
897  trainedNetwork=(TTrainedNetwork*)file2->Get("TTrainedNetwork");
898  cout <<" hid lay 1 size: " << trainedNetwork->getnHiddenLayerSize()[0] << endl;
899  file2->Close();
900  delete file2;
901  */
902 
903  // jn->DumpToFile(name);
904  }
905  }
906 
907  jn->writeNetworkInfo(1);
908  jn->writeNetworkInfo(2);
909  // jn->writeNetworkInfo(3);
910  // jn->writeNetworkInfo(4);
911  // jn->writeNetworkInfo(5);
912 
913 
914  // cout << " Now try to understand how to get the weights..." << endl;
915 
917  Int_t nInput=jn->GetInputDim();
918 
919  cout << " create Trained Network object..." << endl;
920 
921  TTrainedNetwork* trainedNetwork=jn->createTrainedNetwork();
922 
923 /*
924  cout << " now getting value with trained Network ";
925 
926 
927 
928 
929  double inputexample[9]={norm_nVTX(1),
930  norm_nTracksAtVtx(2),
931  norm_nSingleTracks(0),
932  norm_energyFraction(0.6),
933  norm_mass(2500),
934  norm_significance3d(4 ),
935  norm_IP3D(3),
936  norm_cat_pT(3),
937  norm_cat_eta(1)};
938 
939  for (Int_t i=0;i<nInput;++i)
940  {
941  jn->SetInputs(i,inputexample[i]);
942  }
943 
944  cronology.open("weights/trainingCronology.txt",ios_base::app);
945 
946  jn->Evaluate();
947 
948  cronology << "----------------CONSISTENCY CHECK-----------" << endl;
949  cout << "Result 0:" << jn->GetOutput(0);
950  cronology << "Result 0:" << jn->GetOutput(0);
951  cout << " Result 1:" << jn->GetOutput(1);
952  cronology << "Result 0:" << jn->GetOutput(1);
953  cout << " Result 2:" << jn->GetOutput(2) << endl;
954  cronology << " Result 2:" << jn->GetOutput(2) << endl;
955 
956  cout << " Reading back old network " << endl;
957  jn->readBackTrainedNetwork(trainedNetwork);
958 
959  cout <<" resetting input " << endl;
960  for (Int_t i=0;i<nInput;++i)
961  {
962  jn->SetInputs(i,inputexample[i]);
963  }
964 
965  jn->Evaluate();
966 
967  cout << "After reading back - Result 0:" << jn->GetOutput(0);
968  cronology << "After reading back - Result 0:" << jn->GetOutput(0);
969  // << " my: " << result[0] << endl;
970  cout << " After reading back - Result 1:" << jn->GetOutput(1);
971  cronology << "After reading back - Result 1:" << jn->GetOutput(1);
972  //<< " my: " << result[1] << endl;
973  cout << " After reading back - Result 2:" << jn->GetOutput(2) << endl;
974  cronology << "After reading back - Result 2:" << jn->GetOutput(2);
975  // << " my: " << result[2] << endl;
976  */
977 
978  cout << " Now getting histograms from trainingResult" << endl;
979  cronology << " Now getting histograms from trainingResult" << endl;
980 
981  TNetworkToHistoTool myHistoTool;
982 
983  cout << " From network to histo..." << endl;
984  std::vector<TH1*> myHistos=myHistoTool.fromTrainedNetworkToHisto(trainedNetwork);
985 
986  cout << " From histo to network back..." << endl;
987  TTrainedNetwork* trainedNetwork2=myHistoTool.fromHistoToTrainedNetwork(myHistos);
988 
989  cout << " reading back " << endl;
990  jn->readBackTrainedNetwork(trainedNetwork2);
991 
992 /*
993  cout <<" resetting input " << endl;
994  for (Int_t i=0;i<nInput;++i)
995  {
996  jn->SetInputs(i,inputexample[i]);
997  }
998 
999  jn->Evaluate();
1000 
1001  cout << "After reading back - Result 0:" << jn->GetOutput(0);
1002  cronology << "After reading back - Result 0:" << jn->GetOutput(0);
1003  // << " my: " << result[0] << endl;
1004  cout << " After reading back - Result 1:" << jn->GetOutput(1);
1005  cronology << "After reading back - Result 1:" << jn->GetOutput(1);
1006  //<< " my: " << result[1] << endl;
1007  cout << " After reading back - Result 2:" << jn->GetOutput(2) << endl;
1008  cronology << "After reading back - Result 2:" << jn->GetOutput(2);
1009  // << " my: " << result[2] << endl;
1010 
1011 
1012  cout << " Directly from the trainedNetwork read back from HISTOS...!" << endl;
1013 
1014  std::vector<Double_t> inputData;
1015  for (Int_t u=0;u<nInput;++u)
1016  {
1017  inputData.push_back(inputexample[u]);
1018  }
1019 
1020  std::vector<Double_t> outputData=trainedNetwork2->calculateOutputValues(inputData);
1021 
1022  cout << "After reading back - Result 0:" << outputData[0] << endl;
1023  cout << " After reading back - Result 1:" << outputData[1] << endl;
1024  cout << " After reading back - Result 2:" << outputData[2] << endl;
1025 */
1026 
1027 
1028 
1029 
1030  if (epochWithMinimum!=0)
1031  {
1032  cronology << "Minimum stored from Epoch: " << epochWithMinimum << endl;
1033  } else
1034  {
1035  cronology << "Minimum not reached" << endl;
1036  }
1037 
1038  cronology.close();
1039 
1040  if (epochWithMinimum!=0)
1041  {
1042 
1043  TString name=directory;
1044  name+="/Weights";
1045  name+=epochWithMinimum;
1046  name+=".root";
1047  std::cout << " reading back from minimum " << endl;
1048 
1049 
1050  TFile *_file0 = new TFile(name);
1051  TTrainedNetwork* trainedNetwork=(TTrainedNetwork*)_file0->Get("TTrainedNetwork");
1052 
1053  cout << " Reading back network with minimum" << endl;
1054  jn->readBackTrainedNetwork(trainedNetwork);
1055 
1056  TString nameFile=directory;
1057  nameFile+="/weightMinimum.root";
1058 
1059  TFile* file=new TFile(nameFile,"recreate");
1060  trainedNetwork->Write();
1061  file->Write();
1062  file->Close();
1063  delete file;
1064 
1065  cout << " -------------------- " << endl;
1066  cout << " Writing OUTPUT histos " << endl;
1067  TString histoFName=directory;
1068  histoFName+="/histoWeights.root";
1069 
1070  TFile* fileHistos=new TFile(histoFName,"recreate");
1071  TNetworkToHistoTool histoTool;
1072  std::vector<TH1*> myHistos=histoTool.fromTrainedNetworkToHisto(trainedNetwork);
1073  std::vector<TH1*>::const_iterator histoBegin=myHistos.begin();
1074  std::vector<TH1*>::const_iterator histoEnd=myHistos.end();
1075  for (std::vector<TH1*>::const_iterator histoIter=histoBegin;
1076  histoIter!=histoEnd;++histoIter)
1077  {
1078  (*histoIter)->Write();
1079  }
1080  fileHistos->Write();
1081  fileHistos->Close();
1082  delete fileHistos;
1083 
1084  // " filename: " << name << endl;
1085 
1086  // jn->ReadFromFile(name);
1087 
1088  }
1089  else
1090  {
1091  cout << " using network at last iteration (minimum not reached..." << endl;
1092  }
1093 
1094  //here you should create the class... Still open how to deal with this...
1095  // char* myname=const_cast<char*>(static_cast<const char*>(outputclass));
1096  // ierr=mlpsavecf_(myname);
1097 
1098  TString histoTName=directory;
1099  histoTName+="/trainingInfo.root";
1100 
1101  TFile* histoFile=new TFile(histoTName,"recreate");
1102  histoTraining->Write();
1103  histoTesting->Write();
1104  histoFile->Write();
1105  histoFile->Close();
1106  delete histoFile;
1107 
1108  TCanvas* trainingCanvas=new TCanvas("trainingCanvas","trainingCanvas");
1109  histoTraining->SetLineColor(2);
1110  histoTesting->SetLineColor(4);
1111 
1112  histoTraining->GetYaxis()->SetRangeUser(minimumTrain,maximumTrain);
1113  histoTraining->Draw("l");
1114  histoTesting->Draw("lsame");
1115  TString canvasName=directory;
1116  canvasName+="/trainingCurve.eps";
1117  trainingCanvas->SaveAs(canvasName);
1118 
1119 
1120  TCanvas* mlpa_canvas = new TCanvas("jetnet_canvas","Network analysis");
1121  mlpa_canvas->Divide(2,4);
1122 
1123 /*
1124 
1125 // TCanvas* mlpa_canvas_5=gDirectory->Get("mlpa_canvas_5");
1126 // mlpa_canvas_5->SetLogy(kTrue);
1127  gPad->SetLogy();
1128 
1129  // Use the NN to plot the results for each sample
1130  // This will give approx. the same result as DrawNetwork.
1131  // All entries are used, while DrawNetwork focuses on
1132  // the test sample. Also the xaxis range is manually set.
1133  TH1F *bg2 = new TH1F("bg2h", "NN output", 50, -.5, 1.5);
1134  TH1F *bg = new TH1F("bgh", "NN output", 50, -.5, 1.5);
1135  TH1F *sig = new TH1F("sigh", "NN output", 50, -.5, 1.5);
1136 
1137  //sig = 1 part; bg = 2 part; bg2 = 3 part
1138 
1139  TH1F *bg2test = new TH1F("bg2htest", "NN output", 50, -.5, 1.5);
1140  TH1F *bgtest = new TH1F("bghtest", "NN output", 50, -.5, 1.5);
1141  TH1F *sigtest = new TH1F("sightest", "NN output", 50, -.5, 1.5);
1142 
1143  int weight=1;
1144 
1145  for (Int_t i = 0; i < totalN; i++) {
1146 
1147  if (i % 100000 == 0 ) {
1148  std::cout << " First plot. Looping over event " << i << std::endl;
1149  }
1150 
1151  if (i%dilutionFactor!=0&&i%dilutionFactor!=1) continue;
1152 
1153  simu->GetEntry(i);
1154 
1155  for (int u=0;u<sizeX;u++)
1156  {
1157  for (int s=0;s<sizeY;s++)
1158  {
1159  jn->SetInputs( s+u*sizeY, norm_ToT((*matrixOfToT)[u][s]));
1160  }
1161  }
1162  for (int s=0;s<sizeY;s++)
1163  {
1164  jn->SetInputs( sizeX*sizeY+s, norm_pitch((*vectorOfPitchesY)[s]));
1165  }
1166 
1167  jn->SetInputs( (sizeX+1)*sizeY, norm_phi(phi) );
1168  jn->SetInputs( (sizeX+1)*sizeY+1, norm_theta(theta) );
1169 
1170  jn->Evaluate();
1171 
1172  float p1=jn->GetOutput(0);
1173  float p2=jn->GetOutput(1);
1174  float p3=jn->GetOutput(2);
1175 
1176  if (nParticles==1)
1177  {
1178  if (i%dilutionFactor==0)
1179  {
1180  sig->Fill(p1/(p1+p2+p3),weight);
1181  }
1182  else if (i%dilutionFactor==1)
1183  {
1184  sigtest->Fill(p1/(p1+p2+p3),weight);
1185  }
1186  }
1187  if (nParticles==2)
1188  {
1189  if (i%dilutionFactor==0)
1190  {
1191  bg->Fill(p1/(p1+p2+p3),weight);
1192  }
1193  else if (i%dilutionFactor==1)
1194  {
1195  bgtest->Fill(p1/(p1+p2+p3),weight);
1196  }
1197  }
1198  if (nParticles>=3)
1199  {
1200  if (i%dilutionFactor==0)
1201  {
1202  bg2->Fill(p1/(p1+p2+p3),weight);
1203  }
1204  else if (i%dilutionFactor==1)
1205  {
1206  bg2test->Fill(p1/(p1+p2+p3),weight);
1207  }
1208  }
1209  }
1210 
1211  //now you need the maximum
1212  float maximum=1;
1213  for (Int_t a=0;a<bg->GetNbinsX();a++)
1214  {
1215  if (bg->GetBinContent(a)>maximum)
1216  {
1217  maximum=1.2*bg->GetBinContent(a);
1218  }
1219  }
1220 
1221 
1222  bg2->SetLineColor(kYellow);
1223  bg2->SetFillStyle(3008); bg2->SetFillColor(kYellow);
1224  bg->SetLineColor(kBlue);
1225  bg->SetFillStyle(3008); bg->SetFillColor(kBlue);
1226  sig->SetLineColor(kRed);
1227  sig->SetFillStyle(3003); sig->SetFillColor(kRed);
1228  bg2->SetStats(0);
1229  bg->SetStats(0);
1230  sig->SetStats(0);
1231 
1232 
1233  bg2test->SetLineColor(kYellow);
1234  bg2test->SetFillStyle(3008); bg2test->SetFillColor(kYellow);
1235  bgtest->SetLineColor(kBlue);
1236  bgtest->SetFillStyle(3008); bgtest->SetFillColor(kBlue);
1237  sigtest->SetLineColor(kRed);
1238  sigtest->SetFillStyle(3003); sigtest->SetFillColor(kRed);
1239  bg2test->SetStats(0);
1240  bgtest->SetStats(0);
1241  sigtest->SetStats(0);
1242 
1243  mlpa_canvas->cd(1);
1244  gPad->SetLogy();
1245 
1246  bg->GetYaxis()->SetRangeUser(1,maximum);
1247  bgtest->GetYaxis()->SetRangeUser(1,maximum);
1248 
1249  mlpa_canvas->cd(1);
1250  bg->Draw();
1251  bg2->Draw("same");
1252  sig->Draw("same");
1253 
1254  TLegend *legend = new TLegend(.75, .80, .95, .95);
1255  legend->AddEntry(bg2, "particles >=3");
1256  legend->AddEntry(bg, "particles = 2");
1257  legend->AddEntry(sig, "particles = 1");
1258  legend->Draw();
1259 
1260  mlpa_canvas->cd(2);
1261  gPad->SetLogy();
1262 
1263  bgtest->Draw();
1264  bg2test->Draw("same");
1265  sigtest->Draw("same");
1266 
1267  TLegend *legendtest = new TLegend(.75, .80, .95, .95);
1268  legendtest->AddEntry(bg2test, "particles >=3");
1269  legendtest->AddEntry(bgtest, "particles = 2");
1270  legendtest->AddEntry(sigtest, "particles = 1");
1271  legendtest->Draw();
1272 
1273  mlpa_canvas->cd(5);
1274  gPad->SetLogy();
1275  bg->DrawNormalized();
1276  bg2->DrawNormalized("same");
1277  sig->DrawNormalized("same");
1278  legend->Draw();
1279 
1280  mlpa_canvas->cd(6);
1281  gPad->SetLogy();
1282  bgtest->DrawNormalized();
1283  bg2test->DrawNormalized("same");
1284  sigtest->DrawNormalized("same");
1285  legendtest->Draw();
1286 
1287 
1288 
1289  mlpa_canvas->cd(3);
1290  gPad->SetLogy();
1291 
1292  // Use the NN to plot the results for each sample
1293  // This will give approx. the same result as DrawNetwork.
1294  // All entries are used, while DrawNetwork focuses on
1295  // the test sample. Also the xaxis range is manually set.
1296  TH1F *c_bg2 = new TH1F("c_bg2h", "NN output", 50, -.5, 1.5);
1297  TH1F *c_bg = new TH1F("c_bgh", "NN output", 50, -.5, 1.5);
1298  TH1F *c_sig = new TH1F("c_sigh", "NN output", 50, -.5, 1.5);
1299 
1300  TH1F *c_bg2test = new TH1F("c_bg2htest", "NN output", 50, -.5, 1.5);
1301  TH1F *c_bgtest = new TH1F("c_bghtest", "NN output", 50, -.5, 1.5);
1302  TH1F *c_sigtest = new TH1F("c_sightest", "NN output", 50, -.5, 1.5);
1303 
1304  for (Int_t i = 0; i < totalN; i++) {
1305 
1306  if (i % 100000 == 0 ) {
1307  std::cout << " Second plot. Looping over event " << i << std::endl;
1308  }
1309 
1310  if (i%dilutionFactor!=0&&i%dilutionFactor!=1) continue;
1311 
1312  simu->GetEntry(i);
1313 
1314  for (int u=0;u<sizeX;u++)
1315  {
1316  for (int s=0;s<sizeY;s++)
1317  {
1318  jn->SetInputs( s+u*sizeY, norm_ToT((*matrixOfToT)[u][s]));
1319  }
1320  }
1321  for (int s=0;s<sizeY;s++)
1322  {
1323  jn->SetInputs( sizeX*sizeY+s, norm_pitch((*vectorOfPitchesY)[s]));
1324  }
1325 
1326  jn->SetInputs( (sizeX+1)*sizeY, norm_phi(phi) );
1327  jn->SetInputs( (sizeX+1)*sizeY+1, norm_theta(theta) );
1328 
1329  jn->Evaluate();
1330 
1331  float p1=jn->GetOutput(0);
1332  float p2=jn->GetOutput(1);
1333  float p3=jn->GetOutput(2);
1334 
1335  float discr=(p1+p2)/(p1+p2+p3);
1336 
1337  if (nParticles==1)
1338  {
1339  if (i%dilutionFactor==0)
1340  {
1341  c_sig->Fill(discr,weight);
1342  }
1343  else if (i%dilutionFactor==1)
1344  {
1345  c_sigtest->Fill(discr,weight);
1346  }
1347  }
1348  if (nParticles==2)
1349  {
1350  if (i%dilutionFactor==0)
1351  {
1352  c_bg->Fill(discr,weight);
1353  }
1354  else if (i%dilutionFactor==1)
1355  {
1356  c_bgtest->Fill(discr,weight);
1357  }
1358  }
1359  if (nParticles>=3)
1360  {
1361  if (i%dilutionFactor==0)
1362  {
1363  c_bg2->Fill(discr,weight);
1364  }
1365  else if (i%dilutionFactor==1)
1366  {
1367  c_bg2test->Fill(discr,weight);
1368  }
1369  }
1370  }
1371 
1372  //now you need the maximum
1373  maximum=1;
1374  for (Int_t a=0;a<c_bg->GetNbinsX();a++)
1375  {
1376  if (c_bg->GetBinContent(a)>maximum)
1377  {
1378  maximum=1.2*c_bg->GetBinContent(a);
1379  }
1380  }
1381 
1382  c_bg2->SetLineColor(kYellow);
1383  c_bg2->SetFillStyle(3008); c_bg2->SetFillColor(kYellow);
1384  c_bg->SetLineColor(kBlue);
1385  c_bg->SetFillStyle(3008); c_bg->SetFillColor(kBlue);
1386  c_sig->SetLineColor(kRed);
1387  c_sig->SetFillStyle(3003); c_sig->SetFillColor(kRed);
1388  c_bg2->SetStats(0);
1389  c_bg->SetStats(0);
1390  c_sig->SetStats(0);
1391 
1392  c_bg2test->SetLineColor(kYellow);
1393  c_bg2test->SetFillStyle(3008); c_bg2test->SetFillColor(kYellow);
1394  c_bgtest->SetLineColor(kBlue);
1395  c_bgtest->SetFillStyle(3008); c_bgtest->SetFillColor(kBlue);
1396  c_sigtest->SetLineColor(kRed);
1397  c_sigtest->SetFillStyle(3003); c_sigtest->SetFillColor(kRed);
1398  c_bg2test->SetStats(0);
1399  c_bgtest->SetStats(0);
1400  c_sigtest->SetStats(0);
1401 
1402  mlpa_canvas->cd(3);
1403  gPad->SetLogy();
1404 
1405 
1406  c_bg->GetYaxis()->SetRangeUser(1,maximum);
1407  c_bgtest->GetYaxis()->SetRangeUser(1,maximum);
1408 
1409  c_bg->Draw();
1410  c_bg2->Draw("same");
1411  c_sig->Draw("same");
1412 
1413  TLegend *legend2 = new TLegend(.75, .80, .95, .95);
1414  legend2->AddEntry(c_bg2, "particles >=3");
1415  legend2->AddEntry(c_bg, "particles = 2");
1416  legend2->AddEntry(c_sig, "particles = 1");
1417  legend2->Draw();
1418 
1419  mlpa_canvas->cd(4);
1420  gPad->SetLogy();
1421 
1422  c_bgtest->Draw();
1423  c_bg2test->Draw("same");
1424  c_sigtest->Draw("same");
1425 
1426  TLegend *legend2test = new TLegend(.75, .80, .95, .95);
1427  legend2test->AddEntry(c_bg2test, "particles >=3");
1428  legend2test->AddEntry(c_bgtest, "particles = 2");
1429  legend2test->AddEntry(c_sigtest, "particles = 1");
1430  legend2test->Draw();
1431 
1432  mlpa_canvas->cd(7);
1433  gPad->SetLogy();
1434  c_bg->DrawNormalized();
1435  c_bg2->DrawNormalized("same");
1436  c_sig->DrawNormalized("same");
1437  legend2->Draw();
1438 
1439  mlpa_canvas->cd(8);
1440  gPad->SetLogy();
1441  c_bgtest->DrawNormalized();
1442  c_bg2test->DrawNormalized("same");
1443  c_sigtest->DrawNormalized("same");
1444  legend2test->Draw();
1445 
1446 
1447  mlpa_canvas->cd(0);
1448 
1449 
1450  mlpa_canvas->SaveAs("weights/result.eps");
1451 
1452 */
1453 }
1454 
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
trainNN.h
norm_etaModule
double norm_etaModule(const double input)
Definition: NnNormalization.cxx:101
TTrainedNetwork::calculateOutputValues
std::vector< Double_t > calculateOutputValues(std::vector< Double_t > &input) const
Definition: InnerDetector/InDetCalibAlgs/PixelCalibAlgs/NNClusteringCalibration_RunI/TTrainedNetwork.cxx:99
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
TJetNet::GetErrorMeasure
Int_t GetErrorMeasure(void)
Definition: TJetNet.cxx:1156
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
TJetNet::GetLearningRate
Double_t GetLearningRate(void)
Definition: TJetNet.cxx:1177
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
TJetNet::SetPatternsPerUpdate
void SetPatternsPerUpdate(Int_t aValue)
Definition: TJetNet.cxx:1105
DMTest::P
P_v1 P
Definition: P.h:23
extractSporadic.nameFile
string nameFile
Definition: extractSporadic.py:84
TJetNet::readBackTrainedNetwork
void readBackTrainedNetwork(const TTrainedNetwork *)
Definition: TJetNet.cxx:207
TJetNet::writeNetworkInfo
void writeNetworkInfo(Int_t typeOfInfo=0)
Definition: TJetNet.cxx:664
norm_pitch
double norm_pitch(const double input, bool addIBL=false)
Definition: NnNormalization.cxx:32
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:75
TJetNet::GetMomentum
Double_t GetMomentum(void)
Definition: TJetNet.cxx:1182
TJetNet::GetOutput
Double_t GetOutput(Int_t aIndex=0)
Definition: TJetNet.cxx:948
TJetNet::Init
void Init(void)
Definition: TJetNet.cxx:670
python.AthDsoLogger.out
out
Definition: AthDsoLogger.py:71
sigmoid
Double_t sigmoid(Double_t x)
Definition: positions/trainNN.cxx:35
TJetNet::SetOutputTestSet
void SetOutputTestSet(Int_t aPatternInd, Int_t aOutputInd, Double_t aValue)
Definition: TJetNet.h:196
TJetNet::GetLearningRateDecrease
Double_t GetLearningRateDecrease(void)
Definition: TJetNet.cxx:1192
norm_phi
double norm_phi(const double input)
Definition: NnNormalization.cxx:71
TJetNet::SetUpdatingProcedure
void SetUpdatingProcedure(Int_t aValue)
Definition: TJetNet.cxx:1078
badTrackInfo
bool badTrackInfo(bool useTrackEstimate, double theta)
Definition: positions/trainNN.cxx:77
TNetworkToHistoTool
Definition: TNetworkToHistoTool.h:18
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
norm_theta
double norm_theta(const double input)
Definition: NnNormalization.cxx:79
x
#define x
TJetNet::GetUpdatesPerEpoch
Int_t GetUpdatesPerEpoch(void)
Definition: TJetNet.cxx:1146
TJetNet::SetInputTestSet
void SetInputTestSet(Int_t aPatternInd, Int_t aInputInd, Double_t aValue)
Definition: TJetNet.cxx:740
TJetNet::SetEventWeightTrainSet
void SetEventWeightTrainSet(Int_t aPatternInd, Double_t aValue)
Definition: TJetNet.cxx:752
TJetNet::SetMomentum
void SetMomentum(Double_t aValue)
Definition: TJetNet.cxx:1118
TJetNet::Shuffle
void Shuffle(Bool_t aShuffleTrainSet=true, Bool_t aShuffleTestSet=true)
Definition: TJetNet.cxx:1222
norm_posY
double norm_posY(const double input)
Definition: NnNormalization.cxx:118
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
norm_layerNumber
double norm_layerNumber(const double input)
Definition: NnNormalization.cxx:57
TNetworkToHistoTool::fromHistoToTrainedNetwork
TTrainedNetwork * fromHistoToTrainedNetwork(std::vector< TH1 * > &) const
TNetworkToHistoTool::fromTrainedNetworkToHisto
std::vector< TH1 * > fromTrainedNetworkToHisto(TTrainedNetwork *) const
DeMoScan.directory
string directory
Definition: DeMoScan.py:80
lumiFormat.i
int i
Definition: lumiFormat.py:85
TJetNet
Definition: TJetNet.h:41
z
#define z
TJetNet::ReadFromFile
void ReadFromFile(TString aFileName="fort.8")
Definition: TJetNet.cxx:962
norm_ToT
double norm_ToT(const double input)
Definition: NnNormalization.cxx:24
file
TFile * file
Definition: tile_monitor.h:29
TJetNet::SetOutputTrainSet
void SetOutputTrainSet(Int_t aPatternInd, Int_t aOutputInd, Double_t aValue)
Definition: TJetNet.cxx:734
norm_layerType
double norm_layerType(const double input)
Definition: NnNormalization.cxx:64
TJetNet::SetInputTrainSet
void SetInputTrainSet(Int_t aPatternInd, Int_t aInputInd, Double_t aValue)
Definition: TJetNet.cxx:728
TJetNet::Evaluate
void Evaluate(Int_t aPattern)
Definition: TJetNet.cxx:932
TJetNet::GetPatternsPerUpdate
Int_t GetPatternsPerUpdate(void)
Definition: TJetNet.cxx:1172
TTrainedNetwork
Definition: InnerDetector/InDetCalibAlgs/PixelCalibAlgs/NNClusteringCalibration_RunI/TTrainedNetwork.h:21
TJetNet::GetEpochs
Int_t GetEpochs(void)
Definition: TJetNet.h:78
norm_phiBS
double norm_phiBS(const double input)
Definition: NnNormalization.cxx:86
trainNN
void trainNN(TString inputfile, TString outputclass, int nIterations, int dilutionFactor, int nodesFirstLayer, int nodesSecondLayer, int restartTrainingFrom, int nParticlesTraining, bool useTrackEstimate, int nPatternsPerUpdate, double learningRate, double learningRateDecrease, double learningRateMomentum)
Definition: positions/trainNN.cxx:101
TJetNet::GetActivationFunction
Int_t GetActivationFunction(void) const
Definition: TJetNet.cxx:1161
TJetNet::SetErrorMeasure
void SetErrorMeasure(Int_t aValue)
Definition: TJetNet.cxx:1085
TJetNet::SetLearningRate
void SetLearningRate(Double_t aValue)
Definition: TJetNet.cxx:1111
TJetNet::GetOutputDim
Int_t GetOutputDim(void) const
Definition: TJetNet.h:57
TJetNet::SetLearningRateDecrease
void SetLearningRateDecrease(Double_t aValue)
Definition: TJetNet.cxx:1130
skipSingle
bool skipSingle(int nParticles, int iClus, int dilutionFactor)
Definition: positions/trainNN.cxx:56
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
TJetNet::GetInitialWeightsWidth
Double_t GetInitialWeightsWidth(void)
Definition: TJetNet.cxx:1187
TJetNet::SetUpdatesPerEpoch
void SetUpdatesPerEpoch(Int_t aValue)
Definition: TJetNet.cxx:1071
TJetNet::SetInitialWeightsWidth
void SetInitialWeightsWidth(Double_t aValue)
Definition: TJetNet.cxx:1124
norm_posX
double norm_posX(const double input, const bool recenter=false)
Definition: NnNormalization.cxx:108
TJetNet::SetActivationFunction
void SetActivationFunction(Int_t aValue)
Definition: TJetNet.cxx:1091
main
int main()
Definition: positions/trainNN.cxx:86
std::sort
void sort(typename std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, typename std::reverse_iterator< DataModel_detail::iterator< DVL > > end, const Compare &comp)
Specialization of sort for DataVector/List.
Definition: DVL_algorithms.h:623
FullCPAlgorithmsTest_CA.inputfile
dictionary inputfile
Definition: FullCPAlgorithmsTest_CA.py:59
TJetNet::GetInputDim
Int_t GetInputDim(void) const
Definition: TJetNet.h:54
TJetNet::GetHiddenLayerSize
Int_t GetHiddenLayerSize(Int_t number) const
Definition: TJetNet.h:56
TJetNet::createTrainedNetwork
TTrainedNetwork * createTrainedNetwork() const
Definition: TJetNet.cxx:104
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
TJetNet::Test
Double_t Test(void)
Definition: TJetNet.cxx:320
isBadCluster
bool isBadCluster(int sizeX, int nParticles)
Definition: positions/trainNN.cxx:41
test_pyathena.counter
counter
Definition: test_pyathena.py:15
TJetNet::SetEventWeightTestSet
void SetEventWeightTestSet(Int_t aPatternInd, Double_t aValue)
Definition: TJetNet.cxx:758
get_generator_info.command
string command
Definition: get_generator_info.py:38
norm_thetaBS
double norm_thetaBS(const double input)
Definition: NnNormalization.cxx:94
TJetNet::Train
Double_t Train(void)
Definition: TJetNet.cxx:618
fitman.k
k
Definition: fitman.py:528
TJetNet::GetUpdatingProcedure
Int_t GetUpdatingProcedure(void)
Definition: TJetNet.cxx:1151
TJetNet::GetHiddenLayerDim
Int_t GetHiddenLayerDim(void) const
Definition: TJetNet.h:55
Generated on Thu Nov 7 2024 21:29:18 for ATLAS Offline Software by   doxygen 1.8.18