FEDRA emulsion software from the OPERA Collaboration
Public Member Functions | Protected Member Functions | Private Attributes | List of all members
EdbShowAlgID_E_G_Simple Class Reference

#include <EdbShowAlgID_E_G_Simple.h>

Inheritance diagram for EdbShowAlgID_E_G_Simple:
Inheritance graph
[legend]
Collaboration diagram for EdbShowAlgID_E_G_Simple:
Collaboration graph
[legend]

Public Member Functions

 ClassDef (EdbShowAlgID_E_G_Simple, 1)
 
 EdbShowAlgID_E_G_Simple ()
 
 EdbShowAlgID_E_G_Simple (EdbShowerP *shower)
 
 EdbShowAlgID_E_G_Simple (TObjArray *RecoShowerArray)
 
void Execute ()
 
void FindClosestEfficiencyParametrization (Double_t TestAngle, Double_t ReferenceEff)
 
Int_t GetRecoShowerArrayN () const
 
void Print ()
 
void SetCalibrationOffset (Float_t CalibrationOffset)
 
void SetCalibrationSlope (Float_t CalibrationSlope)
 
void SetRecoShowerArray (TObjArray *RecoShowerArray)
 
void SetRecoShowerArrayN (Int_t RecoShowerArrayN)
 
void SetWeightFileString (TString weightstring)
 
virtual ~EdbShowAlgID_E_G_Simple ()
 

Protected Member Functions

void Set0 ()
 

Private Attributes

Float_t eCalibrationOffset
 
Float_t eCalibrationSlope
 
Int_t eParaName
 
TObjArray * eRecoShowerArray
 
Int_t eRecoShowerArrayN
 
TString eWeightFileString
 

Constructor & Destructor Documentation

◆ EdbShowAlgID_E_G_Simple() [1/3]

EdbShowAlgID_E_G_Simple::EdbShowAlgID_E_G_Simple ( )

◆ EdbShowAlgID_E_G_Simple() [2/3]

EdbShowAlgID_E_G_Simple::EdbShowAlgID_E_G_Simple ( EdbShowerP *  shower)
20 {
21  // Default Constructor
22  cout << "EdbShowAlgID_E_G_Simple::EdbShowAlgID_E_G_Simple(EdbShowerP* shower) Default Constructor"<<endl;
23  Set0();
24  eRecoShowerArray=new TObjArray();
25 
26  eRecoShowerArray->Add(shower);
27  SetRecoShowerArrayN(eRecoShowerArray->GetEntries());
28 
29  //cout << "EdbShowAlgID_E_G_Simple::EdbShowAlgID_E_G_Simple(EdbShowerP* shower) Default Constructor Doing Execute right now:"<<endl;
30  //Execute();
31 }
EdbShowAlgID_E_G_Simple::Set0
void Set0()
Definition: EdbShowAlgID_E_G_Simple.cxx:57
EdbShowAlgID_E_G_Simple::SetRecoShowerArrayN
void SetRecoShowerArrayN(Int_t RecoShowerArrayN)
Definition: EdbShowAlgID_E_G_Simple.h:67
EdbShowAlgID_E_G_Simple::eRecoShowerArray
TObjArray * eRecoShowerArray
Definition: EdbShowAlgID_E_G_Simple.h:37

◆ EdbShowAlgID_E_G_Simple() [3/3]

EdbShowAlgID_E_G_Simple::EdbShowAlgID_E_G_Simple ( TObjArray *  RecoShowerArray)
36 {
37  // Default Constructor
38  cout << "EdbShowAlgID_E_G_Simple::EdbShowAlgID_E_G_Simple(TObjArray* RecoShowerArray) Default Constructor"<<endl;
39 
40  Set0();
41  eRecoShowerArray=RecoShowerArray;
42  eRecoShowerArrayN=eRecoShowerArray->GetEntries();
43  cout << "EdbShowAlgID_E_G_Simple::EdbShowAlgID_E_G_Simple(TObjArray* RecoShowerArray) Default Constructor...done."<<endl;
44 }
RecoShowerArray
TObjArray * RecoShowerArray
Definition: Shower_E_FromShowerRoot.C:12
EdbShowAlgID_E_G_Simple::eRecoShowerArrayN
Int_t eRecoShowerArrayN
Definition: EdbShowAlgID_E_G_Simple.h:38

◆ ~EdbShowAlgID_E_G_Simple()

EdbShowAlgID_E_G_Simple::~EdbShowAlgID_E_G_Simple ( )
virtual
49 {
50  // Default Destructor
51  cout << "EdbShowAlgID_E_G_Simple::~EdbShowAlgID_E_G_Simple()"<<endl;
52 }

Member Function Documentation

◆ ClassDef()

EdbShowAlgID_E_G_Simple::ClassDef ( EdbShowAlgID_E_G_Simple  ,
1   
)

◆ Execute()

void EdbShowAlgID_E_G_Simple::Execute ( )

Create ANN Tree and MLP:

70 {
71  cout << "EdbShowAlgID_E_G_Simple::Execute()"<<endl;
72  if (eRecoShowerArrayN==0) {
73  cout << "EdbShowAlgID_E_G_Simple::Execute() eRecoShowerArrayN==0"<<endl;
74  return;
75  }
76 
78  eParaName=2;
79 // TFile* f= new TFile("eee.root","RECREATE");
80  TTree *ANNTree = new TTree("ANNTree", "ANNTree");
81  ANNTree->SetDirectory(0);
82 
83  Float_t inANN[70];
84  ANNTree->Branch("inANN", inANN, "inANN[70]/F");
85 
86  TMultiLayerPerceptron *ANN_MLP;
87  TMultiLayerPerceptron *ANN_MLP_ARRAY[15];
88 
89  TString layout="";
90  int NrOfANNInputNeurons=12;
91 
92  int n_ANN_Neurons_FJ[15];
93  n_ANN_Neurons_FJ[0]=6+10;
94  n_ANN_Neurons_FJ[1]=6+12;
95  n_ANN_Neurons_FJ[2]=6+14;
96  n_ANN_Neurons_FJ[3]=6+16;
97  n_ANN_Neurons_FJ[4]=6+18;
98  n_ANN_Neurons_FJ[5]=6+20;
99  n_ANN_Neurons_FJ[6]=6+23;
100  n_ANN_Neurons_FJ[7]=6+26;
101  n_ANN_Neurons_FJ[8]=6+29;
102  n_ANN_Neurons_FJ[9]=6+32;
103  n_ANN_Neurons_FJ[10]=6+35;
104  n_ANN_Neurons_FJ[11]=6+40;
105  n_ANN_Neurons_FJ[12]=6+45;
106  n_ANN_Neurons_FJ[13]=6+50;
107  n_ANN_Neurons_FJ[14]=6+57;
108 // 10,12,14,16,18,20,23,26,29,32,35,40,45,50
109 
110  for (int k=0; k<15; k++) {
111  layout="";
112  NrOfANNInputNeurons=n_ANN_Neurons_FJ[k];
113  for (Int_t i=1; i<NrOfANNInputNeurons; ++i) layout += "@inANN["+TString(Form("%d",i))+"],";
114  layout += "@inANN["+TString(Form("%d",NrOfANNInputNeurons))+"]:"+TString(Form("%d",NrOfANNInputNeurons+1))+":"+TString(Form("%d",NrOfANNInputNeurons));
115  layout+=":inANN[0]"; // ID CASE
116  //==C== =======================//==C== =======================
117  //cout << "=== DEBUG === NOT ALL ENTRIES FOR ANN TRAINING USED, ONLY A QUICK SOLUTION" << endl;
118  ANN_MLP_ARRAY[k] = new TMultiLayerPerceptron(layout,ANNTree,"","");
119  ANN_MLP=ANN_MLP_ARRAY[k];
120  // if (!quick_run) ANN_MLP = new TMultiLayerPerceptron(layout,ANNInterimTreeInputvariables,"Entry$>5000","Entry$<5000");
121  cout << "=== DEBUG === ANN_MLP->GetStructure() : " << ANN_MLP->GetStructure() << endl;
122  TString ANN_LAYOUT=ANN_MLP_ARRAY[k]->GetStructure();
123  }
124 
125 
126 
127 
128  for (int ll=0; ll<15; ll++) ANN_MLP_ARRAY[ll]->LoadWeights("/home/meisel/Software/ANN_MEGA_ID/weigths_ID_E_GAMMA.txt");
129  cout << "/home/meisel/Software/ANN_MEGA_ID/weigths_ID_E_GAMMA.txt" << endl;
130 
131  // Depending on the Number of Plates for the reconstructed shower,
132  // we choose the right ANN (+layout) and weightfile to be loaded:
133  // Npl() element of [0..10] => layout[0]
134  // Npl() element of [11..12] => layout[1]
135  // Npl() element of [13..14] => layout[2]
136  // Npl() element of [15..16] => layout[3]
137  // Npl() element of [17..18] => layout[4]
138  // Npl() element of [19..20] => layout[5]
139  // Npl() element of [21..23] => layout[6]
140  // Npl() element of [24..26] => layout[7]
141  // Npl() element of [27..29] => layout[8]
142  // Npl() element of [30..32] => layout[9]
143  // Npl() element of [33..35] => layout[10]
144  // Npl() element of [36..40] => layout[11]
145  // Npl() element of [41..45] => layout[12]
146  // Npl() element of [46..50] => layout[13]
147 
148  TH1F* hist_E = new TH1F("hist_E","hist_E",100,0,10);
149 
150  for (int i=0; i<eRecoShowerArrayN; i++) {
151 // for (int i=0; i<4; i++) {
152 
153  EdbShowerP* shower=(EdbShowerP*) eRecoShowerArray->At(i);
154  shower->Print();
155 
156  shower->BuildParametrisation_FJ();
157  EdbShowerP::Para_FJ eShowAlgE_Para_FJ=shower->GetPara_FJ();
158 // Para_FJ eShowAlgE_Para_FJ=shower->GetPara_FJ();
159  shower->PrintParametrisation_FJ();
160 
161  cout << "shower->Npl() " << shower->Npl() << endl;
162  int check_Npl= shower->Npl();
163  int check_Npl_index=0;
164  if (check_Npl<=10) check_Npl_index=0;
165  if (check_Npl>=11 && check_Npl<=12) check_Npl_index=1;
166  if (check_Npl>=13 && check_Npl<=14) check_Npl_index=2;
167  if (check_Npl>=15 && check_Npl<=16) check_Npl_index=3;
168  if (check_Npl>=17 && check_Npl<=18) check_Npl_index=4;
169  if (check_Npl>=19 && check_Npl<=20) check_Npl_index=5;
170  if (check_Npl>=21 && check_Npl<=23) check_Npl_index=6;
171  if (check_Npl>=24 && check_Npl<=26) check_Npl_index=7;
172  if (check_Npl>=27 && check_Npl<=29) check_Npl_index=8;
173  if (check_Npl>=30 && check_Npl<=32) check_Npl_index=9;
174  if (check_Npl>=33 && check_Npl<=35) check_Npl_index=10;
175  cout << "check_Npl_index " << check_Npl_index << endl;
176 
177 
178 // cout << " =========================== FORCE check_Npl_index TO BE 20 PLATES !!! ======================="<<endl;
179 // check_Npl_index=5;
180 // cout << " =========================== check_Npl_index=5; !!! ======================="<<endl;
181 
182 // cout << " =========================== FORCE check_Npl_index TO BE 32 PLATES !!! ======================="<<endl;
183 // check_Npl_index=9;
184 // cout << " =========================== check_Npl_index=9; !!! ======================="<<endl;
185 
186 
187 
188  if (check_Npl>20) check_Npl_index=9;
189  if (check_Npl<=20) check_Npl_index=5;
190  if (check_Npl<=10) check_Npl_index=0;
191 
192 
193  check_Npl_index=5;
194 
195  cout << " =========================== check_Npl_index= " << check_Npl_index << " ; !!! ======================="<<endl;
196 
197  NrOfANNInputNeurons=n_ANN_Neurons_FJ[check_Npl_index];
198  ANN_MLP=ANN_MLP_ARRAY[check_Npl_index];
199  cout << "Using the following layout: " << ANN_MLP->GetStructure() << endl;
200  ANN_MLP->Print();
201 
202  // inANN[0] is ALWAYS reserved for the quantity value to be estimated
203  // (E,Id,...)
204 
205  inANN[1]=eShowAlgE_Para_FJ.ShowerAxisAngle;
206  inANN[2]=eShowAlgE_Para_FJ.nseg;
207  inANN[3]=eShowAlgE_Para_FJ.BT_deltaR_mean;
208  inANN[4]=eShowAlgE_Para_FJ.BT_deltaR_rms;
209  inANN[5]=eShowAlgE_Para_FJ.BT_deltaT_mean;
210  inANN[6]=eShowAlgE_Para_FJ.BT_deltaT_rms;
211  for (int ii=0; ii<57; ii++) {
212  cout << setw(3) << eShowAlgE_Para_FJ.longprofile[ii];
213  inANN[7+ii]= eShowAlgE_Para_FJ.longprofile[ii];
214  }
215  cout << endl;
216  // Fill Tree:
217  ANNTree->Fill();
218 
219  Double_t params[70];
220  Double_t val;
221  Double_t EnergyCalibration_Offset=0.8;
222  Double_t EnergyCalibration_Slope=0.94;
223 
224 
225  EnergyCalibration_Offset=0.0;
226  EnergyCalibration_Slope=1.0;
227 
228 
229  for (Int_t j=1; j<=NrOfANNInputNeurons; ++j) {
230  params[j-1]=inANN[j];
231  cout << "=== DEBUG === : j params[j-1]=inANN[j] " << j<< " " << params[j-1] << endl;
232  }
233  val=(ANN_MLP->Evaluate(0,params));
234  cout << "DEBUG: i, val (=mlp1->Evaluate(0,params)) ======================= real energy " << i << " " << val << " == " << inANN[0] << endl;
235 
236  if ( ::isnan(val) ) cout << " val is NAN " << endl;
237 
238 
239  val=(val-EnergyCalibration_Offset)/EnergyCalibration_Slope;
240  cout << "---- With Correctiono factor val=(val-EnergyCalibration_Offset)/EnergyCalibration_Slope: " << EnergyCalibration_Offset << " " << EnergyCalibration_Slope << " " << val<< endl;
241 
242 
243 
244  // Quick Estimation of the sigma (from ANN_MEGA_ENERGY)
245  Double_t sigma=0.83*TMath::Sqrt(val)-0.45;
246 
247  cout << "Estimated Energy = " << val << " +- " << sigma << "..." << endl;
248 
249 
250 // FindClosestEfficiencyParametrization(inANN[1],0.85);
251 
252 // hist_E->Fill(val);
253  }
254 
255 // hist_E->Draw();
256 
257  return;
258 }
EdbShowAlgID_E_G_Simple::eParaName
Int_t eParaName
Definition: EdbShowAlgID_E_G_Simple.h:39
EdbShowerP
Definition: EdbShowerP.h:28
EdbShowerP::Print
void Print()
-— EXPERIMENTAL STATUS —
Definition: EdbShowerP.cxx:2333
EdbShowerP::Npl
Int_t Npl() const
Definition: EdbShowerP.h:409
EdbShowerP::BuildParametrisation_FJ
void BuildParametrisation_FJ()
Definition: EdbShowerP.cxx:1235
EdbShowerP::PrintParametrisation_FJ
void PrintParametrisation_FJ()
Definition: EdbShowerP.cxx:2392
EdbShowerP::GetPara_FJ
Para_FJ GetPara_FJ() const
Definition: EdbShowerP.h:551
EdbShowerP::Para_FJ
Definition: EdbShowerP.h:143
EdbShowerP::Para_FJ::BT_deltaT_rms
Float_t BT_deltaT_rms
Definition: EdbShowerP.h:150
EdbShowerP::Para_FJ::nseg
Int_t nseg
Definition: EdbShowerP.h:146
EdbShowerP::Para_FJ::BT_deltaR_rms
Float_t BT_deltaR_rms
Definition: EdbShowerP.h:148
EdbShowerP::Para_FJ::BT_deltaR_mean
Float_t BT_deltaR_mean
Definition: EdbShowerP.h:147
EdbShowerP::Para_FJ::longprofile
Int_t longprofile[57]
Definition: EdbShowerP.h:151
EdbShowerP::Para_FJ::ShowerAxisAngle
Float_t ShowerAxisAngle
Definition: EdbShowerP.h:145
EdbShowerP::Para_FJ::BT_deltaT_mean
Float_t BT_deltaT_mean
Definition: EdbShowerP.h:149
LoadWeights
mlp LoadWeights("testweights.txt")
params
Double_t params[3]
Definition: testBGReduction_By_ANN.C:84

◆ FindClosestEfficiencyParametrization()

void EdbShowAlgID_E_G_Simple::FindClosestEfficiencyParametrization ( Double_t  TestAngle = 0.0,
Double_t  ReferenceEff = 1.0 
)
264 {
265 
266  cout << "TestAngle " << TestAngle << endl;
267  cout << "ReferenceEff " << ReferenceEff << endl;
268 
269  TF1 *all = new TF1("all","1.0-0.00000001*x",0,0.95);
270  TF1 *edefault = new TF1("default","0.94-0.216*x-0.767*x*x+1.865*x*x*x",0,0.95);
271  TF1 *elletroni = new TF1("elletroni","0.79+0.38*x-7.63*x*x+25.13*x*x*x-24.6*x*x*x*x",0,0.95);
272  TF1 *neuchmicro = new TF1("neuchmicro","0.94-0.955*x+1.80*x*x-0.95*x*x*x",0,0.95);
273  TF1 *MiddleFix = new TF1("MiddleFix","0.5*(0.888361-1.299*x+1.49198*x*x+1.64668*x*x*x-2.63758*x*x*x*x+0.79+0.38*x-7.63*x*x+25.13*x*x*x-24.6*x*x*x*x)",0,0.95);
274  TF1 *LowEff = new TF1("LowEff","0.85*0.5*(0.888361-1.299*x+1.49198*x*x+1.64668*x*x*x-2.63758*x*x*x*x+0.79+0.38*x-7.63*x*x+25.13*x*x*x-24.6*x*x*x*x)",0,0.95);
275 
276 // Double_t TestAngle=0.05;
277 // Double_t ReferenceEff=0.776;
278 
279  cout << "all->Eval(TestAngle) "<< all->Eval(TestAngle)<< endl;
280  cout << "LowEff->Eval(TestAngle) "<< LowEff->Eval(TestAngle)<< endl;
281  cout << "MiddleFix->Eval(TestAngle) "<< MiddleFix->Eval(TestAngle)<< endl;
282  cout << "elletroni->Eval(TestAngle) "<< elletroni->Eval(TestAngle)<< endl;
283  cout << "neuchmicro->Eval(TestAngle) "<< neuchmicro->Eval(TestAngle)<< endl;
284  cout << "default->Eval(TestAngle) "<< edefault->Eval(TestAngle)<< endl;
285 
286 //Measure distance of angle to estimated angle:
287 
288  Double_t dist[9];
289 
290  dist[0]=TMath::Abs(ReferenceEff-edefault->Eval(TestAngle));
291  dist[1]=TMath::Abs(ReferenceEff-neuchmicro->Eval(TestAngle));
292  dist[2]=TMath::Abs(ReferenceEff-elletroni->Eval(TestAngle));
293  dist[3]=TMath::Abs(ReferenceEff-MiddleFix->Eval(TestAngle));
294  dist[4]=TMath::Abs(ReferenceEff-LowEff->Eval(TestAngle));
295  dist[5]=TMath::Abs(ReferenceEff-all->Eval(TestAngle));
296 
297  Double_t dist_min=1;
298  Int_t best_i=-1;
299  for (int i=0; i<6; i++) if (abs(dist[i])<dist_min) {
300  dist_min=dist[i];
301  best_i=i;
302  }
303 
304  cout << "Miminum distance = " << dist_min << endl;
305  cout << "Best numerb efficiency = ";
306  if (best_i==0) cout << " default " <<endl;
307  if (best_i==1) cout << " neuchmicro " <<endl;
308  if (best_i==2) cout << " elletroni " <<endl;
309  if (best_i==3) cout << " MiddleFix " <<endl;
310  if (best_i==4) cout << " LowEff " <<endl;
311  if (best_i==5) cout << " all " <<endl;
312 
313  return;
314 
315 }

◆ GetRecoShowerArrayN()

Int_t EdbShowAlgID_E_G_Simple::GetRecoShowerArrayN ( ) const
inline
71  {
72  return eRecoShowerArrayN;
73  }

◆ Print()

void EdbShowAlgID_E_G_Simple::Print ( )

◆ Set0()

void EdbShowAlgID_E_G_Simple::Set0 ( )
protected
58 {
59  eCalibrationOffset=0;
60  eCalibrationSlope=1;
61  eRecoShowerArray=NULL;
62  eRecoShowerArrayN=0;
63  eWeightFileString="NULL";
64  return;
65 }
EdbShowAlgID_E_G_Simple::eCalibrationOffset
Float_t eCalibrationOffset
Definition: EdbShowAlgID_E_G_Simple.h:42
EdbShowAlgID_E_G_Simple::eCalibrationSlope
Float_t eCalibrationSlope
Definition: EdbShowAlgID_E_G_Simple.h:43
EdbShowAlgID_E_G_Simple::eWeightFileString
TString eWeightFileString
Definition: EdbShowAlgID_E_G_Simple.h:46
NULL
#define NULL
Definition: nidaqmx.h:84

◆ SetCalibrationOffset()

void EdbShowAlgID_E_G_Simple::SetCalibrationOffset ( Float_t  CalibrationOffset)
inline
76  {
77  eCalibrationOffset=CalibrationOffset;
78  }

◆ SetCalibrationSlope()

void EdbShowAlgID_E_G_Simple::SetCalibrationSlope ( Float_t  CalibrationSlope)
inline
79  {
80  eCalibrationSlope=CalibrationSlope;
81  }

◆ SetRecoShowerArray()

void EdbShowAlgID_E_G_Simple::SetRecoShowerArray ( TObjArray *  RecoShowerArray)
inline
63  {
64  eRecoShowerArray=RecoShowerArray;
65  eRecoShowerArrayN=eRecoShowerArray->GetEntriesFast();
66  }

◆ SetRecoShowerArrayN()

void EdbShowAlgID_E_G_Simple::SetRecoShowerArrayN ( Int_t  RecoShowerArrayN)
inline
67  {
68  eRecoShowerArrayN=RecoShowerArrayN;
69  }

◆ SetWeightFileString()

void EdbShowAlgID_E_G_Simple::SetWeightFileString ( TString  weightstring)
320 {
321  eWeightFileString=weightstring;
322  cout << "EdbShowAlgID_E_G_Simple::SetWeightFileString() eWeightFileString = " << eWeightFileString << endl;
323  return;
324 }

Member Data Documentation

◆ eCalibrationOffset

Float_t EdbShowAlgID_E_G_Simple::eCalibrationOffset
private

◆ eCalibrationSlope

Float_t EdbShowAlgID_E_G_Simple::eCalibrationSlope
private

◆ eParaName

Int_t EdbShowAlgID_E_G_Simple::eParaName
private

◆ eRecoShowerArray

TObjArray* EdbShowAlgID_E_G_Simple::eRecoShowerArray
private

◆ eRecoShowerArrayN

Int_t EdbShowAlgID_E_G_Simple::eRecoShowerArrayN
private

◆ eWeightFileString

TString EdbShowAlgID_E_G_Simple::eWeightFileString
private
The documentation for this class was generated from the following files: