CbmRoot
CbmRichElectronIdAnn.cxx
Go to the documentation of this file.
1 
8 #include "CbmRichElectronIdAnn.h"
9 #include "CbmGlobalTrack.h"
10 #include "CbmRichGeoManager.h"
11 #include "CbmRichRing.h"
12 #include "CbmRichUtil.h"
13 #include "FairLogger.h"
14 #include "FairRootManager.h"
15 #include "TClonesArray.h"
16 #include "TMath.h"
17 #include "TMultiLayerPerceptron.h"
18 #include "TSystem.h"
19 #include "TTree.h"
20 
21 #include <iostream>
22 
23 using std::cout;
24 using std::endl;
25 
27  : fAnnWeights(""), fNN(NULL), fGlobalTracks(NULL), fRichRings(NULL) {
28  Init();
29 }
30 
32 
34  if (fNN != NULL) { delete fNN; }
35 
36  if (CbmRichGeoManager::GetInstance().fGP->fGeometryType
38  fAnnWeights = string(gSystem->Getenv("VMCWORKDIR"))
39  + "/parameters/rich/rich_v17a_elid_ann_weights.txt";
40  } else if (CbmRichGeoManager::GetInstance().fGP->fGeometryType
42  fAnnWeights = string(gSystem->Getenv("VMCWORKDIR"))
43  + "/parameters/rich/rich_v16a_elid_ann_weights.txt";
44  } else {
45  fAnnWeights = string(gSystem->Getenv("VMCWORKDIR"))
46  + "/parameters/rich/rich_v17a_elid_ann_weights.txt";
47  }
48 
49  TTree* simu = new TTree("MonteCarlo", "MontecarloData");
50  Double_t x[9];
51  Double_t xOut;
52 
53  simu->Branch("x0", &x[0], "x0/D");
54  simu->Branch("x1", &x[1], "x1/D");
55  simu->Branch("x2", &x[2], "x2/D");
56  simu->Branch("x3", &x[3], "x3/D");
57  simu->Branch("x4", &x[4], "x4/D");
58  simu->Branch("x5", &x[5], "x5/D");
59  simu->Branch("x6", &x[6], "x6/D");
60  simu->Branch("x7", &x[7], "x7/D");
61  simu->Branch("x8", &x[8], "x8/D");
62  simu->Branch("xOut", &xOut, "xOut/D");
63 
64  fNN = new TMultiLayerPerceptron("x0,x1,x2,x3,x4,x5,x6,x7,x8:18:xOut", simu);
65  cout << "-I- CbmRichElIdAnn: get NeuralNet weight parameters from: "
66  << fAnnWeights << endl;
67  fNN->LoadWeights(fAnnWeights.c_str());
68 
69  FairRootManager* ioman = FairRootManager::Instance();
70  if (ioman != NULL) {
71  fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
72  if (fRichRings == NULL) {
73  LOG(error) << "CbmRichElectronIdAnn::Init() fRichRings == NULL";
74  }
75  fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
76  if (fGlobalTracks == NULL) {
77  LOG(error) << "CbmRichElectronIdAnn::Init() fGlobalTracks == NULL";
78  }
79  } else {
80  LOG(error) << "FairRootManager::Instance() == NULL";
81  }
82 }
83 
84 double CbmRichElectronIdAnn::CalculateAnnValue(int globalTrackIndex,
85  double momentum) {
86  double errorValue = -1.;
87  if (globalTrackIndex < 0) return errorValue;
88 
89  if (fGlobalTracks == NULL || fRichRings == NULL) return -1;
90 
91  const CbmGlobalTrack* globalTrack =
92  static_cast<const CbmGlobalTrack*>(fGlobalTracks->At(globalTrackIndex));
93  if (globalTrack == NULL) return errorValue;
94 
95  Int_t richId = globalTrack->GetRichRingIndex();
96  if (richId == -1) return errorValue;
97  const CbmRichRing* richRing =
98  static_cast<const CbmRichRing*>(fRichRings->At(richId));
99  if (richRing == NULL) return errorValue;
100 
101  double rtDistance = CbmRichUtil::GetRingTrackDistance(globalTrackIndex);
102 
103  if (richRing->GetAaxis() >= 10. || richRing->GetAaxis() <= 0.
104  || richRing->GetBaxis() >= 10. || richRing->GetBaxis() <= 0.
105  || richRing->GetNofHits() <= 5. || rtDistance <= 0. || rtDistance >= 999.
106  || richRing->GetRadialPosition() <= 0.
107  || richRing->GetRadialPosition() >= 999. || richRing->GetPhi() <= -6.5
108  || richRing->GetPhi() >= 6.5 || richRing->GetRadialAngle() <= -6.5
109  || richRing->GetRadialAngle() >= 6.5) {
110 
111  return -1.;
112  }
113  double params[9];
114  params[0] = richRing->GetAaxis() / 10.;
115  params[1] = richRing->GetBaxis() / 10.;
116  params[2] = (richRing->GetPhi() + 1.57) / 3.14;
117  params[3] = richRing->GetRadialAngle() / 6.28;
118  params[4] = (richRing->GetChi2() / richRing->GetNDF()) / 1.2;
119  params[5] = richRing->GetRadialPosition() / 110.;
120  params[6] = richRing->GetNofHits() / 45.;
121  params[7] = rtDistance / 5.;
122  params[8] = momentum / 15.;
123 
124  for (int k = 0; k < 9; k++) {
125  if (params[k] < 0.0) params[k] = 0.0;
126  if (params[k] > 1.0) params[k] = 1.0;
127  }
128 
129  double nnEval = fNN->Evaluate(0, params);
130 
131  if (TMath::IsNaN(nnEval) == 1) {
132  //cout << " -W- CbmRichElectronIdAnn: nnEval nan " << endl;
133  nnEval = -1.;
134  }
135 
136  return nnEval;
137 }
CbmRichGeometryTypeTwoWings
@ CbmRichGeometryTypeTwoWings
Definition: CbmRichRecGeoPar.h:24
CbmRichGeoManager::GetInstance
static CbmRichGeoManager & GetInstance()
Definition: CbmRichGeoManager.h:29
CbmRichElectronIdAnn::fRichRings
TClonesArray * fRichRings
Definition: CbmRichElectronIdAnn.h:73
CbmGlobalTrack::GetRichRingIndex
Int_t GetRichRingIndex() const
Definition: CbmGlobalTrack.h:41
CbmRichRing::GetChi2
Double_t GetChi2() const
Definition: CbmRichRing.h:95
CbmGlobalTrack.h
CbmRichRing::GetNofHits
Int_t GetNofHits() const
Definition: CbmRichRing.h:40
CbmRichRing
Definition: CbmRichRing.h:17
CbmRichRing.h
CbmRichRing::GetRadialPosition
Float_t GetRadialPosition() const
Definition: CbmRichRing.cxx:181
CbmRichGeoManager.h
CbmRichUtil::GetRingTrackDistance
static Double_t GetRingTrackDistance(Int_t globalTrackId)
Definition: alignment/CbmRichUtil.h:20
CbmRichRing::GetRadialAngle
Double_t GetRadialAngle() const
Definition: CbmRichRing.cxx:189
CbmRichElectronIdAnn::CbmRichElectronIdAnn
CbmRichElectronIdAnn()
Standard constructor.
Definition: CbmRichElectronIdAnn.cxx:26
CbmRichGeometryTypeCylindrical
@ CbmRichGeometryTypeCylindrical
Definition: CbmRichRecGeoPar.h:25
CbmRichElectronIdAnn::fAnnWeights
string fAnnWeights
Set path to the file with ANN weights.
Definition: CbmRichElectronIdAnn.h:69
CbmRichRing::GetAaxis
Double_t GetAaxis() const
Definition: CbmRichRing.h:83
CbmRichRing::GetBaxis
Double_t GetBaxis() const
Definition: CbmRichRing.h:84
CbmGlobalTrack
Definition: CbmGlobalTrack.h:26
CbmRichRing::GetPhi
Double_t GetPhi() const
Definition: CbmRichRing.h:87
CbmRichElectronIdAnn::fNN
TMultiLayerPerceptron * fNN
Definition: CbmRichElectronIdAnn.h:70
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmRichElectronIdAnn::CalculateAnnValue
double CalculateAnnValue(int globalTrackIndex, double momentum)
Calculate output value of the ANN.
Definition: CbmRichElectronIdAnn.cxx:84
CbmRichRing::GetNDF
Double_t GetNDF() const
Definition: CbmRichRing.h:96
CbmRichElectronIdAnn::fGlobalTracks
TClonesArray * fGlobalTracks
Definition: CbmRichElectronIdAnn.h:72
CbmRichElectronIdAnn.h
Implementation of the electron identification algorithm in the RICH detector using Artificial Neural ...
CbmRichElectronIdAnn::~CbmRichElectronIdAnn
virtual ~CbmRichElectronIdAnn()
Destructor.
Definition: CbmRichElectronIdAnn.cxx:31
CbmRichElectronIdAnn::Init
void Init()
Initialize ANN before use.
Definition: CbmRichElectronIdAnn.cxx:33