ATLAS Offline Software
validateNNnum.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 <TH1F.h>
9 #include <TLegend.h>
10 #include <iostream>
11 #include <TPad.h>
12 //#include <stdio.h>
13 #include <string.h>
14 //#include <stdlib.h>
15 #include <math.h>
16 #include "../TJetNet.h"
17 #include "../doNormalization.C"
18 #include "Riostream.h"
19 #include "../TNetworkToHistoTool.h"
20 #include <TSystem.h>
21 
22 #include "../TTrainedNetwork.h"
23 
24 #include <sstream>
25 #include "TMatrixD.h"
26 #include "TVectorD.h"
27 //#include "trainNN.h"
28 #include <algorithm>
29 #include <vector>
30 #include <utility>
31 #include <TLatex.h>
32 #include <TChain.h>
33 #include <TGraphErrors.h>
34 
35 using namespace std;
36 
37 Double_t sigmoid(Double_t x)
38 {
39  return 1./(1.+exp(-2*x));
40 }
41 
42 using namespace std;
43 
44 
45 bool isBadCluster( int sizeX, int nParticles ){
46 
47 
48  if(sizeX==-100){return true;}
49  if(nParticles==0){return true;}
50  if(nParticles!=1 && nParticles!=2 && nParticles!=3){return true;}
51 
52 
53 
54  return false;
55 
56 
57 }
58 
59 
60 double refineRMS(TH1D* hist ){
61 
62  double HistoRms = hist->GetRMS();
63  double HistMean = hist->GetMean();
64  //hist->SetAxisRange(HistMean - 3.*HistoRms, HistMean + 3.*HistoRms );
65  double RMS = hist->GetRMS();
66 
67  return RMS;
68 
69 }
70 
71 
72 
73 
74 void validateNN( bool useTrackEstimate ) {
75 
76 
77 
78  // gROOT->SetStyle("Plain");
79 
80  TString inputfile="root://castoratlas//castor/cern.ch/user/l/lorenzi/TrkValidation17.0.0.5_J5.root";
81  // TString inputfile="root://castoratlas//castor/cern.ch/user/l/lorenzi/TrkValidation17.0.0.5_J5_onTrack.root";
82  // TString inputfile="/tmp/lorenzi/user.lorenzi.J5.TRKVAL.NNtrain.onTrack1.110620162506/TrkValidation_J5_onTrack.root";
83  // TString inputfile="/tmp/lorenzi/user.lorenzi.T1.TRKVAL.NNtrain.onTrack.110620153532/TrkValidation_T1_onTrack.root";
84 
85  // TString inputfile="/tmp/lorenzi/17.0.0.5/TrkValidation.root";
86 
87  TFile *file= TFile::Open(inputfile);
88 
89  TTree *simu = (TTree*)file->Get("Validation/NNinput");
90 
91 
92 
93  std::cout << " Training sample obtained... " << std::endl;
94 
95  vector<int> *NN_sizeX;
96  vector<int> *NN_sizeY;
97  vector<vector<float> > *NN_matrixOfToT;
98  vector<vector<float> > *NN_vectorOfPitchesY;
99  vector<int> *NN_ClusterPixLayer;
100  vector<int> *NN_ClusterPixBarrelEC;
101  vector<float> *NN_phiBS;
102  vector<float> *NN_thetaBS;
103  vector<float> *NN_etaModule;
104  vector<bool> *NN_useTrackInfo;
105  vector<int> *NN_columnWeightedPosition;
106  vector<int> *NN_rowWeightedPosition;
107  vector<double> *NN_localColumnWeightedPosition;
108  vector<double> *NN_localRowWeightedPosition;
109 
110  vector<vector<float> > *NN_positionX;
111  vector<vector<float> > *NN_positionY;
112  vector<vector<float> > *NN_position_idX;
113  vector<vector<float> > *NN_position_idY;
114  vector<vector<float> > *NN_theta;
115  vector<vector<float> > *NN_phi;
116 
117  // List of branches
118  TBranch *b_NN_sizeX;
119  TBranch *b_NN_sizeY;
120  TBranch *b_NN_matrixOfToT;
121  TBranch *b_NN_vectorOfPitchesY;
122  TBranch *b_NN_ClusterPixLayer;
123  TBranch *b_NN_ClusterPixBarrelEC;
124  TBranch *b_NN_phiBS;
125  TBranch *b_NN_thetaBS;
126  TBranch *b_NN_etaModule;
127  TBranch *b_NN_useTrackInfo;
128  TBranch *b_NN_columnWeightedPosition;
129  TBranch *b_NN_rowWeightedPosition;
130  TBranch *b_NN_localColumnWeightedPosition;
131  TBranch *b_NN_localRowWeightedPosition;
132  TBranch *b_NN_positionX;
133  TBranch *b_NN_positionY;
134  TBranch *b_NN_position_idX;
135  TBranch *b_NN_position_idY;
136  TBranch *b_NN_theta;
137  TBranch *b_NN_phi;
138 
139 
140 
141  NN_sizeX = 0;
142  NN_sizeY = 0;
143  NN_matrixOfToT = 0;
144  NN_vectorOfPitchesY = 0;
145  NN_ClusterPixLayer = 0;
146  NN_ClusterPixBarrelEC = 0;
147  NN_phiBS = 0;
148  NN_thetaBS = 0;
149  NN_etaModule = 0;
150  NN_useTrackInfo = 0;
151  NN_columnWeightedPosition = 0;
152  NN_rowWeightedPosition = 0;
153  NN_localColumnWeightedPosition = 0;
154  NN_localRowWeightedPosition = 0;
155  NN_positionX = 0;
156  NN_positionY = 0;
157  NN_position_idX = 0;
158  NN_position_idY = 0;
159  NN_theta = 0;
160  NN_phi = 0;
161  // Set branch addresses and branch pointers
162  // if (!tree) return 0;
163  // TTree* simu = tree;
164  // fCurrent = -1;
165  simu->SetMakeClass(1);
166 
167  simu->SetBranchAddress("NN_sizeX", &NN_sizeX, &b_NN_sizeX);
168  simu->SetBranchAddress("NN_sizeY", &NN_sizeY, &b_NN_sizeY);
169  simu->SetBranchAddress("NN_matrixOfToT", &NN_matrixOfToT, &b_NN_matrixOfToT);
170  simu->SetBranchAddress("NN_vectorOfPitchesY", &NN_vectorOfPitchesY, &b_NN_vectorOfPitchesY);
171  simu->SetBranchAddress("NN_ClusterPixLayer", &NN_ClusterPixLayer, &b_NN_ClusterPixLayer);
172  simu->SetBranchAddress("NN_ClusterPixBarrelEC", &NN_ClusterPixBarrelEC, &b_NN_ClusterPixBarrelEC);
173  simu->SetBranchAddress("NN_phiBS", &NN_phiBS, &b_NN_phiBS);
174  simu->SetBranchAddress("NN_thetaBS", &NN_thetaBS, &b_NN_thetaBS);
175  simu->SetBranchAddress("NN_etaModule", &NN_etaModule, &b_NN_etaModule);
176  simu->SetBranchAddress("NN_useTrackInfo", &NN_useTrackInfo, &b_NN_useTrackInfo);
177  simu->SetBranchAddress("NN_columnWeightedPosition", &NN_columnWeightedPosition, &b_NN_columnWeightedPosition);
178  simu->SetBranchAddress("NN_rowWeightedPosition", &NN_rowWeightedPosition, &b_NN_rowWeightedPosition);
179 
180  simu->SetBranchAddress("NN_localColumnWeightedPosition", &NN_localColumnWeightedPosition, &b_NN_localColumnWeightedPosition);
181  simu->SetBranchAddress("NN_localRowWeightedPosition", &NN_localRowWeightedPosition, &b_NN_localRowWeightedPosition);
182 
183  simu->SetBranchAddress("NN_positionX", &NN_positionX, &b_NN_positionX);
184  simu->SetBranchAddress("NN_positionY", &NN_positionY, &b_NN_positionY);
185  simu->SetBranchAddress("NN_position_idX", &NN_position_idX, &b_NN_position_idX);
186  simu->SetBranchAddress("NN_position_idY", &NN_position_idY, &b_NN_position_idY);
187 
188  simu->SetBranchAddress("NN_theta", &NN_theta, &b_NN_theta);
189  simu->SetBranchAddress("NN_phi", &NN_phi, &b_NN_phi);
190 
191 
192 
193 
194  cout << "Branches set..." << endl;
195 
196 // int nParticlesTraining = 1 ;
197 
198 
199  TString pathWithoutTracks = "/afs/cern.ch/user/g/giacinto/scratch0/PixelClusterisationTF/jetNet/jetNetJune2011/withoutTracks/Weights.root";
200  TString pathWithTracks = "/afs/cern.ch/user/g/giacinto/scratch0/PixelClusterisationTF/prepareJetNetNtupleApril2011/withTracks/Weights.root";
201 
202 
203  TString name;
204 
205  if(!useTrackEstimate){
206 
207  name+= pathWithoutTracks;
208 
209  }else{
210 
211  name+=pathWithTracks;
212 
213  }
214 
215 
216  // getting a postion trained network from file
217  TFile *_file0 = new TFile(name);
218  TTrainedNetwork* numberTrainedNetwork=(TTrainedNetwork*)_file0->Get("TTrainedNetwork");
219 
220  cout << " Reading back network with minimum" << endl;
221 
222 
223  Int_t sizeX=-7;
224  Int_t sizeY=-7;
225 
226 
227  int iClus=0;
228 
229  TH1D* histo1 = new TH1D("1","1",1000,0,1);
230  TH1D* histo2 = new TH1D("2","2",1000,0,1);
231  TH1D* histo3 = new TH1D("3","3",1000,0,1);
232 
233  TH1D* ahisto2= new TH1D("a2","a2",1000,0,1);
234  TH1D* ahisto3= new TH1D("a3","a3",1000,0,1);
235 
236  // Loop over entries:
237  for (Int_t i = 0; i < simu->GetEntries(); i++) {
238 
239  if (i % 1000 == 0 ) {
240  std::cout << " Counting training / testing events in sample. Looping over event " << i << " cluster: " << iClus << std::endl;
241  }
242 
243  if (i > 50000 ) break;
244  // if(iClus > 10000) break;
245 
246  simu->GetEntry(i);
247 
248  for( unsigned int clus =0; clus<NN_sizeX->size(); clus++ ){
249  // cout << clus << endl;
250  vector<float> *matrixOfToT=0;
251  vector<float> *vectorOfPitchesY=0;
252 
253  Float_t phiBS;
254  Float_t thetaBS;
255  Float_t etaModule;
256  Int_t ClusterPixLayer;
257  Int_t ClusterPixBarrelEC;
258 
259  std::vector<float> * positionX=0;
260  std::vector<float> * positionY=0;
261 
262  std::vector<float> * position_idX=0;
263  std::vector<float> * position_idY=0;
264 
265 
266  std::vector<float> * thetaTr=0;
267  std::vector<float> * phiTr=0;
268 
269  double localColumnWeightedPosition;// = 0;
270  double localRowWeightedPosition;// = 0;
271 
272  double columnWeightedPosition;// = 0;
273  double rowWeightedPosition;// = 0;
274 
275 
276 
277  sizeX = (*NN_sizeX)[clus];
278  sizeY = (*NN_sizeY)[clus];
279 
280  matrixOfToT=&(*NN_matrixOfToT)[clus];
281  vectorOfPitchesY=&(*NN_vectorOfPitchesY)[clus];
282 
283  phiBS = (*NN_phiBS)[clus];
284  thetaBS =(*NN_thetaBS)[clus];
285  etaModule =(*NN_etaModule)[clus];
286 
287  ClusterPixLayer=(*NN_ClusterPixLayer)[clus];
288  ClusterPixBarrelEC = (*NN_ClusterPixBarrelEC)[clus];
289 
290  positionX =&(*NN_positionX)[clus];
291  positionY =&(*NN_positionY)[clus];
292 
293  position_idX =&(*NN_position_idX)[clus];
294  position_idY =&(*NN_position_idY)[clus];
295 
296  thetaTr = &(*NN_theta)[clus];
297  phiTr = &(*NN_phi)[clus];
298 
299 
300  localColumnWeightedPosition =(*NN_localColumnWeightedPosition)[clus];
301  localRowWeightedPosition =(*NN_localRowWeightedPosition)[clus];
302 
303  columnWeightedPosition =(*NN_columnWeightedPosition)[clus];
304  rowWeightedPosition =(*NN_rowWeightedPosition)[clus];
305 
306  // if( ClusterPixLayer!=0 || ClusterPixBarrelEC!=0) continue;
307 
308 
309 
310  unsigned int nParticles = positionX->size();
311 
312 
313 
314  if (isBadCluster(sizeX, nParticles ) )continue;
315 
316  thetaTr = &(*NN_theta)[clus];
317  phiTr = &(*NN_phi)[clus];
318 
319 
320 
321  // loop over the particles;
322  for( unsigned int P = 0; P < positionX->size(); P++){
323  iClus++;
324 
325  if(iClus > 5000) break;
326 
327  double theta = (*thetaTr)[P];
328  double phi = (*phiTr)[P];
329 
330  std::vector<Double_t> inputData;
331 
332  for(unsigned int ME =0; ME < matrixOfToT->size(); ME++){
333 
334  inputData.push_back(norm_ToT((*matrixOfToT)[ME]));
335  }
336 
337  for (int s=0;s<sizeY;s++)
338  {
339  inputData.push_back(norm_pitch((*vectorOfPitchesY)[s]));
340  }
341 
342  inputData.push_back(norm_layerNumber(ClusterPixLayer));
343  inputData.push_back(norm_layerType(ClusterPixBarrelEC));
344 
345  if (useTrackEstimate)
346  {
347  inputData.push_back(norm_phi(phi));
348  inputData.push_back(norm_theta(theta));
349  }
350  else
351  {
352  inputData.push_back(norm_phiBS(phiBS));
353  inputData.push_back(norm_thetaBS(thetaBS));
354  inputData.push_back(norm_etaModule(etaModule));
355  }
356  // cout << "formatted" << endl;
357 
358 
359  vector<double> outputNN = numberTrainedNetwork->calculateOutputValues(inputData);
360 
361  int nParticles= positionX->size();
362 
363  if( nParticles ==1 ) {
364 
365  histo1->Fill(outputNN[0]/(outputNN[0]+outputNN[1]+outputNN[2]));
366 
367 
368  }
369 
370  if( nParticles ==2 ) {
371 
372  histo2->Fill(outputNN[0]/(outputNN[0]+outputNN[1]+outputNN[2]));
373 
374  if( outputNN[0]/(outputNN[0]+outputNN[1]+outputNN[2])<0.98 ){
375 
376  ahisto2->Fill( outputNN[1]/(outputNN[1]+outputNN[2]) );
377  }
378 
379  }
380 
381  if( nParticles ==3 ) {
382  histo3->Fill(outputNN[0]/(outputNN[0]+outputNN[1]+outputNN[2]));
383 
384  if( outputNN[0]/(outputNN[0]+outputNN[1]+outputNN[2])<0.98 ) {
385 
386  ahisto3->Fill(outputNN[1]/(outputNN[1]+outputNN[2]) );
387 
388  }
389 
390  }
391 
392 
393 
394 
395 
396  }// end loop over th particles
397  }// end loop over cluster
398  }// end Loop over entries
399 
400 
401 
402 
403 
404 
405  double* integral1= histo1->GetIntegral();
406  double* integral2= histo2->GetIntegral();
407  double* integral3= histo3->GetIntegral();
408 
409  double xcut[1000];
410 
411  double integral4[1000];
412  double integral5[1000];
413  for (int u=0;u<1000;u++)
414  {
415  xcut[u]=histo2->GetBinCenter(u);
416  integral4[u]=1-integral2[u];
417  integral5[u]=1-integral3[u];
418  }
419 
420  TGraphErrors* Have1Vs2=new TGraphErrors(1000,integral1,integral4,0,0);
421  TGraphErrors* Have1Vs3=new TGraphErrors(1000,integral1,integral5,0,0);
422  TGraphErrors* HaveCutVs2=new TGraphErrors(1000,xcut,integral4,0,0);
423 
424  TCanvas * c1=new TCanvas("c1","c1",600,600);
425 
426  Have1Vs2->SetLineColor(2);
427  Have1Vs2->SetMarkerColor(2);
428  Have1Vs2->SetMarkerSize(0.4);
429  Have1Vs2->SetMarkerStyle(3);
430  Have1Vs2->Draw("AP");
431 
432  Have1Vs2->GetXaxis()->SetTitle("Fraction of wrongly split single-cluster");
433  Have1Vs2->GetYaxis()->SetTitle("Fraction of not split multi-clusters");
434 
435 
436 
437 
438  double* aintegral2= ahisto2->GetIntegral();
439  double* aintegral3= ahisto3->GetIntegral();
440 
441  double axcut[1000];
442 
443  double aintegral4[1000];
444  double aintegral5[1000];
445  for (int u=0;u<1000;u++)
446  {
447  axcut[u]=ahisto2->GetBinCenter(u);
448  aintegral4[u]=1-aintegral2[u];
449  aintegral5[u]=1-aintegral3[u];
450  }
451 
452 // double* aintegral21= ahisto21->GetIntegral();
453 
454 
455  TCanvas* c4=new TCanvas("c4","c4",600,600);
456 
457  TGraphErrors* aHave1Vs3=new TGraphErrors(1000,aintegral2,aintegral5,0,0);
458  TGraphErrors* aHaveCutVs2=new TGraphErrors(1000,axcut,aintegral4,0,0);
459 
460 
461  aHave1Vs3->GetXaxis()->SetTitle("Fraction of wrongly split cluster (2 #rightarrow 3)");
462  aHave1Vs3->GetYaxis()->SetTitle("Fraction of not split multi-clusters (3 #rightarrow 2)");
463 
464 
465  aHave1Vs3->SetLineColor(4);
466  aHave1Vs3->SetMarkerColor(4);
467  aHave1Vs3->SetMarkerSize(0.5);
468  aHave1Vs3->SetMarkerStyle(3);
469  aHave1Vs3->Draw("AP");
470 
471 
472  // TCanvas * c2=new TCanvas("c2","c2",600,600);
473  // HaveCutVs2->Draw("AP");
474 
475 
476 
477 
478 
479 
480 
481 }
482 
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
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
plotmaker.hist
hist
Definition: plotmaker.py:148
norm_pitch
double norm_pitch(const double input, bool addIBL=false)
Definition: NnNormalization.cxx:32
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:71
extractSporadic.c1
c1
Definition: extractSporadic.py:134
norm_phi
double norm_phi(const double input)
Definition: NnNormalization.cxx:71
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
norm_theta
double norm_theta(const double input)
Definition: NnNormalization.cxx:79
x
#define x
Trk::u
@ u
Enums for curvilinear frames.
Definition: ParamDefs.h:83
norm_layerNumber
double norm_layerNumber(const double input)
Definition: NnNormalization.cxx:57
sigmoid
Double_t sigmoid(Double_t x)
Definition: validateNNnum.cxx:37
lumiFormat.i
int i
Definition: lumiFormat.py:92
norm_ToT
double norm_ToT(const double input)
Definition: NnNormalization.cxx:24
file
TFile * file
Definition: tile_monitor.h:29
norm_layerType
double norm_layerType(const double input)
Definition: NnNormalization.cxx:64
TTrainedNetwork
Definition: InnerDetector/InDetCalibAlgs/PixelCalibAlgs/NNClusteringCalibration_RunI/TTrainedNetwork.h:21
norm_phiBS
double norm_phiBS(const double input)
Definition: NnNormalization.cxx:86
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
refineRMS
double refineRMS(TH1D *hist)
Definition: validateNNnum.cxx:60
FullCPAlgorithmsTest_CA.inputfile
dictionary inputfile
Definition: FullCPAlgorithmsTest_CA.py:59
DiTauMassTools::HistInfoV2::RMS
@ RMS
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:31
isBadCluster
bool isBadCluster(int sizeX, int nParticles)
Definition: validateNNnum.cxx:45
validateNN
void validateNN(bool useTrackEstimate)
Definition: validateNNnum.cxx:74
norm_thetaBS
double norm_thetaBS(const double input)
Definition: NnNormalization.cxx:94
Generated on Thu Jul 4 2024 21:32:22 for ATLAS Offline Software by   doxygen 1.8.18