25 #include "TClonesArray.h"
26 #include "TMultiLayerPerceptron.h"
28 #include <boost/assign/list_of.hpp>
42 using boost::assign::list_of;
51 fRichRingMatches(NULL),
53 fStsTrackMatches(NULL),
57 fMinNofHitsInRichRing(7),
59 fMaxNofTrainSamples(5000),
101 for (
int i = 0; i < 2; i++){
102 if (i == 0) ss =
"Electron";
103 if (i == 1) ss =
"Pion";
105 fhAaxis[i] =
new TH1D(
string(
"fhAaxis"+ss).c_str(),
"fhAaxis;A axis [cm];Counter", 50, 0., 8.);
107 fhBaxis[i] =
new TH1D(
string(
"fhBAxis"+ss).c_str(),
"fhBAxis;B axis [cm];Counter", 50, 0., 8.);
113 fhDistTrueMatch[i] =
new TH1D(
string(
"fhDistTrueMatch"+ss).c_str(),
"fhDistTrueMatch;Ring-track distance [cm];Counter", 50, 0., 5.);
115 fhDistMisMatch[i] =
new TH1D(
string(
"fhDistMisMatch"+ss).c_str(),
"fhDistMisMatch;Ring-track distance [cm];Counter", 50, 0., 5.);
117 fhNofHits[i] =
new TH1D(
string(
"fhNofHits"+ss).c_str(),
"fhNofHits;Number of hits;Counter", 50, 0, 50);
119 fhChi2[i] =
new TH1D(
string(
"fhChi2"+ss).c_str(),
"fhChi2;#Chi^{2};Counter", 100, 0., 1.);
121 fhRadPos[i] =
new TH1D(
string(
"fhRadPos"+ss).c_str(),
"fhRadPos;Radial position [cm];Counter", 150, 0., 150.);
123 fhAaxisVsMom[i] =
new TH2D(
string(
"fhAaxisVsMom"+ss).c_str(),
"fhAaxisVsMom;P [GeV/c];A axis [cm]",30, 0, 15, 50, 0, 10);
125 fhBaxisVsMom[i] =
new TH2D(
string(
"fhBAxisVsMom"+ss).c_str(),
"fhBAxisVsMom;P [GeV/c];B axis [cm]",30, 0, 15, 50, 0, 10);
127 fhPhiVsRadAng[i] =
new TH2D(
string(
"fhPhiVsRadAng"+ss).c_str(),
"fhPhiVsRadAng;Phi [rad];Radial angle [rad]", 50, -2. ,2.,50, 0. , 6.3);
130 fhAnnOutput[i] =
new TH1D(
string(
"fhAnnOutput"+ss).c_str(),
"ANN output;ANN output;Counter",100, -1.2, 1.2);
132 fhCumProb[i] =
new TH1D(
string(
"fhCumProb"+ss).c_str(),
"ANN output;ANN output;Cumulative probability",100, -1.2, 1.2);
144 cout <<
"InitStatus CbmRichTrainAnnElectrons::Init()"<<endl;
147 if (NULL == ioman) { Fatal(
"CbmRichTrainAnnElectrons::Init",
"RootManager not instantised!");}
150 if ( NULL ==
fRichHits) { Fatal(
"CbmRichTrainAnnElectrons::Init",
"No RichHit array!");}
153 if ( NULL ==
fRichRings) { Fatal(
"CbmRichTrainAnnElectrons::Init",
"No RichRing array!");}
156 if ( NULL ==
fRichPoints) { Fatal(
"CbmRichTrainAnnElectrons::Init",
"No RichPoint array!");}
159 if ( NULL ==
fMCTracks) { Fatal(
"CbmRichTrainAnnElectrons::Init",
"No MCTrack array!");}
162 if ( NULL ==
fRichRingMatches) { Fatal(
"CbmRichTrainAnnElectrons::Init",
"No RichRingMatch array!");}
165 if ( NULL ==
fRichProj) { Fatal(
"CbmRichTrainAnnElectrons::Init",
"No RichProjection array!");}
168 if ( NULL ==
fStsTrackMatches) { Fatal(
"CbmRichTrainAnnElectrons::Init",
"No track match array!");}
171 if ( NULL ==
fGlobalTracks) { Fatal(
"CbmRichTrainAnnElectrons::Init",
"No global track array!");}
174 if ( NULL ==
fStsTracks) { Fatal(
"CbmRichTrainAnnElectrons::Init",
"No STSTrack array!");}
182 cout << endl <<
"-I- CbmRichTrainAnnElectrons, event " <<
fEventNum << endl;
185 cout <<
"Nof Electrons = " <<
fRElIdParams[0].size() << endl;
186 cout <<
"Nof Pions = " <<
fRElIdParams[1].size() << endl;
193 for (Int_t iTrack=0; iTrack < nGlTracks; iTrack++) {
195 if (NULL == gTrack)
continue;
197 if (stsIndex == -1)
continue;
198 CbmStsTrack* stsTrack = (CbmStsTrack*)
fStsTracks->At(stsIndex);
199 if (NULL == stsTrack)
continue;
201 if (NULL == stsTrackMatch)
continue;
205 if (richIndex == -1)
continue;
207 if (NULL == ring)
continue;
209 if (NULL == richRingMatch)
continue;
213 if (NULL == track)
continue;
223 Double_t lPercTrue = 0;
224 if (lFoundHits >= 3){
225 lPercTrue = (Double_t)richRingMatch->
GetNofTrueHits() / (Double_t)lFoundHits;
227 Bool_t isTrueFound = (lPercTrue >
fQuota);
241 if (pdg == 11 && motherId == -1 && isTrueFound &&
242 mcIdSts == mcIdRich && mcIdRich != -1){
257 if (pdg == 11 && motherId == -1 && isTrueFound &&
258 mcIdSts != mcIdRich && mcIdRich != -1){
264 if ( pdg == 211 && mcIdRich != -1){
269 if (mcIdRich == mcIdSts) {
288 TTree *simu =
new TTree (
"MonteCarlo",
"MontecarloData");
292 simu->Branch(
"x0", &x[0],
"x0/D");
293 simu->Branch(
"x1", &x[1],
"x1/D");
294 simu->Branch(
"x2", &x[2],
"x2/D");
295 simu->Branch(
"x3", &x[3],
"x3/D");
296 simu->Branch(
"x4", &x[4],
"x4/D");
297 simu->Branch(
"x5", &x[5],
"x5/D");
298 simu->Branch(
"x6", &x[6],
"x6/D");
299 simu->Branch(
"x7", &x[7],
"x7/D");
300 simu->Branch(
"x8", &x[8],
"x8/D");
301 simu->Branch(
"xOut", &xOut,
"xOut/D");
303 for (
int j = 0; j < 2; j++){
315 for (
int k = 0;
k < 9;
k++){
316 if (x[
k] < 0.0) x[
k] = 0.0;
317 if (x[
k] > 1.0) x[
k] = 1.0;
320 if (j == 0) xOut = 1.;
321 if (j == 1) xOut = -1.;
327 TMultiLayerPerceptron network(
"x0,x1,x2,x3,x4,x5,x6,x7,x8:18:xOut",simu,
"Entry$+1");
329 network.Train(300,
"text,update=10");
330 network.DumpWeights(
"rich_elid_ann_weights.txt");
338 for (
int j = 0; j < 2; j++){
350 for (
int k = 0;
k < 9;
k++){
351 if (params[
k] < 0.0) params[
k] = 0.0;
352 if (params[
k] > 1.0) params[
k] = 1.0;
355 Double_t netEval = network.Evaluate(0,params);
369 cout <<
"nof electrons = " <<
fRElIdParams[0].size() << endl;
370 cout <<
"nof pions = " <<
fRElIdParams[1].size() << endl;
373 cout <<
"ANN cut = " <<
fAnnCut << endl;
375 Double_t cumProbFake = 0.;
376 Double_t cumProbTrue = 0.;
377 Int_t nofFake = (Int_t)
fhAnnOutput[1]->GetEntries();
378 Int_t nofTrue = (Int_t)
fhAnnOutput[0]->GetEntries();
379 for (Int_t i = 1; i <=
fhAnnOutput[1]->GetNbinsX(); i++){
381 fhCumProb[1]->SetBinContent(i, (Double_t)cumProbFake/nofFake);
384 fhCumProb[0]->SetBinContent(i, 1. - (Double_t)cumProbTrue / nofTrue);
388 TCanvas*
c1 =
new TCanvas(
"ann_electrons_ann_output",
"ann_electrons_ann_output", 500, 500);
392 TCanvas*
c2 =
new TCanvas(
"ann_electrons_cum_prob",
"ann_electrons_cum_prob", 500, 500);
397 TCanvas*
c3 =
new TCanvas(
"ann_electrons_params_ab",
"ann_electrons_params_ab", 1200, 600);
409 TCanvas* c3_2 =
new TCanvas(
"ann_electrons_params_1",
"ann_electrons_params_1", 1500, 600);
424 list_of(
"e^{#pm} true match")(
"e^{#pm} mis match")(
"#pi^{#pm} true match")(
"#pi^{#pm} mis match"),
428 TCanvas* c3_1 =
new TCanvas(
"ann_electrons_params_2",
"ann_electrons_params_2", 1500, 600);
438 TCanvas*
c4 =
new TCanvas(
"ann_electrons_params_2d",
"ann_electrons_params_2d", 600, 900);
468 for (
int i = 0; i <
fHists.size(); i++){