diff --git a/OptimizationMLw.C b/OptimizationMLw.C
new file mode 100644
index 0000000000000000000000000000000000000000..3558c3945099c5376549b0d2eab59b1fd76f7fcf
--- /dev/null
+++ b/OptimizationMLw.C
@@ -0,0 +1,536 @@
+/* Copyright (C) 2019-2020 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Anna Senger [committer] */
+
+//---------------------------------------------------
+//
+// invariant mass spectra for omega mesons @ 8-10 A GeV/c
+//
+// CbmAnaDimuonAnalysis has to be used:
+// - with UseCuts(kFALSE) for type 0 (cuts optimization)
+// - with UseCuts(kTRUE) for type 1 (geometry optimization)
+//
+// Invariant mass spectrum for background is calculated using super-event technics
+//
+// Anna Senger a.senger@gsi.de
+//
+//--------------------------------------------------
+// ------ [ for ML part ] ----
+// Authors: Pawan Kumar Sharma, Raktim Mukerjee, Abhishek Sharma, Partha Pratim Bhaduri
+// Selection of Muon track candidates using Machine learnig classifier from TMVA
+//---------------------------------------------------
+
+#include "TMVA/Tools.h" //
+#include "TMVA/Reader.h" //
+#include "TMVA/MethodCuts.h" //
+#include <map>
+
+using namespace TMVA; //
+
+void OptimizationMLw(Int_t energy = 8, Int_t NofFiles = 100, Int_t type = 0, TString version = "v23a",
+ TString dir = "../../../data/v23a/sis100_muon_lmvm")
+
+{
+ gStyle->SetCanvasColor(10);
+ gStyle->SetFrameFillColor(10);
+ gStyle->SetHistLineWidth(4);
+ gStyle->SetPadColor(10);
+ gStyle->SetStatColor(10);
+ gStyle->SetPalette(55);
+ gStyle->SetPaintTextFormat("2.1f");
+
+ //----------- This loads the library from TMVA ----------------------
+ TMVA::Tools::Instance();
+ TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );
+
+ Float_t v1,v2,v3,v4,v5,v6,v7,v8,v9;
+
+ reader->AddVariable("NofMUCH",&v1);
+ reader->AddVariable("NofSTS",&v2);
+ reader->AddVariable("NofTRD",&v3);
+ reader->AddVariable("Chi2MUCH",&v4);
+ reader->AddVariable("Chi2STS",&v5);
+ reader->AddVariable("Chi2Vertex",&v6);
+ reader->AddVariable("Chi2TRD",&v7);
+ reader->AddVariable("M",&v8);
+ reader->AddVariable("P",&v9);
+
+ // Book the MVA methods
+
+ TString MLdir = "../../../../v23aML/part/Data_v23a_omega/weights/";
+ TString prefix = "TMVAClassification";
+ TString methodName="BDTG-1";
+ TString weightfile=MLdir+prefix+"_BDTG-1.weights.xml";
+ reader->BookMVA(methodName, weightfile);
+
+
+ // Double_t sgn_mult = 2.28721 * 7.4 * 1e-5; // for 8AGeV from CBM-wiki (thermal model)
+ Double_t sgn_mult = 19 * 9 * 1e-5; // for 10AGeV from HSD code
+
+ TH1D* invM_sgn = new TH1D("invM_omega", "invM_omega", 400, 0, 4);
+
+ (invM_sgn->GetXaxis())->SetTitle("m_{inv} (GeV/c^{2})");
+ (invM_sgn->GetYaxis())->SetTitle("counts/(events #times 10 MeV/c^{2})");
+ invM_sgn->SetLineColor(kRed);
+ invM_sgn->SetMarkerColor(kRed);
+ invM_sgn->SetMarkerStyle(20);
+
+ TH1D* invM_bg = (TH1D*) invM_sgn->Clone("invM_bg");
+ invM_bg->SetTitle("invM_bg");
+ invM_bg->SetLineColor(kBlack);
+ invM_bg->SetMarkerColor(kBlack);
+
+ TH1D* invM_full = (TH1D*) invM_sgn->Clone("invM_full");
+ invM_full->SetTitle("invM_full");
+ invM_full->SetLineColor(kBlack);
+ invM_full->SetMarkerColor(kBlack);
+
+
+ TClonesArray* MuPlus = new TClonesArray("CbmAnaMuonCandidate");
+ TClonesArray* MuMinus = new TClonesArray("CbmAnaMuonCandidate");
+ TTree* InputTree;
+
+ Double_t NofEvents[2];
+
+ Int_t NEvents_sgn=0, NEvents_bg=0;
+
+ //-------------- soft cuts ---------------------------------
+
+ Int_t MUCHhits = 8;
+ Int_t STShits = 4;
+ Int_t TRDhits = 1;
+ Int_t TOFhits = 1;
+
+ Double_t MUCHchi2 = 20.0;
+ Double_t STSchi2 = 20.0;
+ Double_t Vchi2 = 20;
+
+ Int_t sigmaTOFCut = 2;
+
+ TString dir_par = getenv("VMCWORKDIR");
+ TString name = dir_par + "/parameters/much/TOF8gev_fitParam_sigma2.sis100_muon_lmvm.root";
+
+ Double_t p0min, p1min, p2min;
+ Double_t p0max, p1max, p2max;
+
+ TFile* FF = new TFile(name);
+ TTree* MinParam = (TTree*) FF->Get("MinParam");
+ MinParam->SetBranchAddress("p0", &p0min);
+ MinParam->SetBranchAddress("p1", &p1min);
+ MinParam->SetBranchAddress("p2", &p2min);
+ MinParam->GetEntry(0);
+
+ TTree* MaxParam = (TTree*) FF->Get("MaxParam");
+ MaxParam->SetBranchAddress("p0", &p0max);
+ MaxParam->SetBranchAddress("p1", &p1max);
+ MaxParam->SetBranchAddress("p2", &p2max);
+ MaxParam->GetEntry(0);
+
+ //--------------------------------------------------------
+
+ for (int k = 1; k < NofFiles + 1; k++) {
+ name.Form("%s/%dgev/omega/%d/muons.ana.root", dir.Data(), energy, k);
+
+ TFile* f = new TFile(name);
+ if (f->IsZombie() || f->GetNkeys() < 1 || f->TestBit(TFile::kRecovered)) {
+ f->Close();
+ continue;
+ }
+
+ if (k % 100 == 0) cout << "Input File " << k << " : " << f->GetName() << endl;
+
+ InputTree = (TTree*) f->Get("cbmsim");
+ if (!InputTree) {
+ f->Close();
+ continue;
+ }
+
+ InputTree->SetBranchAddress("MuPlus", &MuPlus);
+ InputTree->SetBranchAddress("MuMinus", &MuMinus);
+
+ //-----------------------------------------------------------
+
+ for (int iEvent = 0; iEvent < InputTree->GetEntries(); iEvent++) {
+ InputTree->GetEntry(iEvent);
+
+ NofEvents[0]++;
+ NEvents_sgn++;
+
+ int NofPlus = MuPlus->GetEntriesFast();
+ int NofMinus = MuMinus->GetEntriesFast();
+
+ if (NofPlus == 0 || NofMinus == 0) continue;
+
+ for (int iPart = 0; iPart < NofPlus; iPart++) {
+ CbmAnaMuonCandidate* mu_pl = (CbmAnaMuonCandidate*) MuPlus->At(iPart);
+
+ if (type == 0) {
+
+
+ if (mu_pl->GetNStsHits() < STShits) continue;
+ if (mu_pl->GetNMuchHits() < MUCHhits) continue;
+ if (mu_pl->GetNTrdHits() < TRDhits) continue;
+ if (mu_pl->GetNTofHits() < TOFhits) continue;
+
+ if (mu_pl->GetChiToVertex() > Vchi2) continue;
+ if (mu_pl->GetChiSts() > STSchi2) continue;
+ if (mu_pl->GetChiMuch() > MUCHchi2) continue;
+
+ v1 = mu_pl->GetNMuchHits();
+ v2 = mu_pl->GetNStsHits();
+ v3 = mu_pl->GetNTrdHits();
+ v4 = mu_pl->GetChiMuch();
+ v5 = mu_pl->GetChiSts();
+ v6 = mu_pl->GetChiToVertex();
+ v7 = mu_pl->GetChiTrd();
+ v8 = mu_pl->GetTofM();
+ v9 = mu_pl->GetMomentum()->P();
+
+ double y1 = reader->EvaluateMVA("BDTG-1"); // ML evaluation
+ if (y1 < .72) continue; // Using a cut on the ML output
+
+ }
+
+ TLorentzVector* P_pl = mu_pl->GetMomentum();
+
+ Double_t mass = mu_pl->GetTofM();
+ Double_t mom = P_pl->P();
+
+
+ if (type == 0)
+ if (mass < (p0min + p1min * mom + p2min * mom * mom) || mass > (p0max + p1max * mom + p2max * mom * mom)) continue;
+
+
+
+ for (int jPart = 0; jPart < NofMinus; jPart++) {
+ CbmAnaMuonCandidate* mu_mn = (CbmAnaMuonCandidate*) MuMinus->At(jPart);
+
+
+ if (type == 0) {
+
+
+ if (mu_mn->GetNStsHits() < STShits) continue;
+ if (mu_mn->GetNMuchHits() < MUCHhits) continue;
+ if (mu_mn->GetNTrdHits() < TRDhits) continue;
+ if (mu_mn->GetNTofHits() < TOFhits) continue;
+
+ if (mu_mn->GetChiToVertex() > Vchi2) continue;
+ if (mu_mn->GetChiSts() > STSchi2) continue;
+ if (mu_mn->GetChiMuch() > MUCHchi2) continue;
+
+
+ v1 = mu_mn->GetNMuchHits();
+ v2 = mu_mn->GetNStsHits();
+ v3 = mu_mn->GetNTrdHits();
+ v4 = mu_mn->GetChiMuch();
+ v5 = mu_mn->GetChiSts();
+ v6 = mu_mn->GetChiToVertex();
+ v7 = mu_mn->GetChiTrd();
+ v8 = mu_mn->GetTofM();
+ v9 = mu_mn->GetMomentum()->P();
+
+ double y1 = reader->EvaluateMVA( "BDTG-1"); // ML evaluation
+ if (y1 < .72) continue; // Using a cut on the ML output
+
+ }
+
+ TLorentzVector* P_mn = mu_mn->GetMomentum();
+
+ mass = mu_mn->GetTofM();
+ mom = P_mn->P();
+
+ if (type == 0)
+ if (mass < (p0min + p1min * mom + p2min * mom * mom) || mass > (p0max + p1max * mom + p2max * mom * mom)) continue;
+
+ TLorentzVector M(*P_pl + *P_mn);
+ invM_sgn->Fill(M.M(), sgn_mult / 10.); // normalized with bin value
+
+
+ }
+ }
+ }
+ f->Close();
+ }
+
+ TLorentzVector P1, P2, M;
+ TTree* Plus = new TTree("Plus", "part-");
+ Plus->Branch("P0", &P1(0), "Px/D:Py:Pz:E", 10000000);
+
+ TTree* Minus = new TTree("Minus", "part+");
+ Minus->Branch("P0", &P2(0), "Px/D:Py:Pz:E", 10000000);
+
+ for (int k = 1; k < NofFiles + 1; k++) {
+ name.Form("%s/%dgev/centr/%d/muons.ana.root", dir.Data(), energy, k);
+
+ TFile* f = new TFile(name);
+ if (f->IsZombie() || f->GetNkeys() < 1 || f->TestBit(TFile::kRecovered)) {
+ f->Close();
+ continue;
+ }
+
+ if (k % 100 == 0) cout << "Input File " << k << " : " << f->GetName() << endl;
+
+ InputTree = (TTree*) f->Get("cbmsim");
+ if (!InputTree) {
+ f->Close();
+ continue;
+ }
+
+ InputTree->SetBranchAddress("MuPlus", &MuPlus);
+ InputTree->SetBranchAddress("MuMinus", &MuMinus);
+
+ //-----------------------------------------------------------
+
+ for (int iEvent = 0; iEvent < InputTree->GetEntries(); iEvent++) {
+ InputTree->GetEntry(iEvent);
+ // std::cout<<" background analysis started "<<std::endl;
+ NofEvents[1]++;
+ NEvents_bg++;
+
+ int NofPlus = MuPlus->GetEntriesFast();
+ int NofMinus = MuMinus->GetEntriesFast();
+
+ for (int iPart = 0; iPart < NofPlus; iPart++) {
+ CbmAnaMuonCandidate* mu_pl = (CbmAnaMuonCandidate*) MuPlus->At(iPart);
+ if (type == 0) {
+ if (mu_pl->GetNStsHits() < STShits) continue;
+ if (mu_pl->GetNMuchHits() < MUCHhits) continue;
+ if (mu_pl->GetNTrdHits() < TRDhits) continue;
+ if (mu_pl->GetNTofHits() < TOFhits) continue;
+
+ if (mu_pl->GetChiToVertex() > Vchi2) continue;
+ if (mu_pl->GetChiSts() > STSchi2) continue;
+ if (mu_pl->GetChiMuch() > MUCHchi2) continue;
+
+ v1 = mu_pl->GetNMuchHits();
+ v2 = mu_pl->GetNStsHits();
+ v3 = mu_pl->GetNTrdHits();
+ v4 = mu_pl->GetChiMuch();
+ v5 = mu_pl->GetChiSts();
+ v6 = mu_pl->GetChiToVertex();
+ v7 = mu_pl->GetChiTrd();
+ v8 = mu_pl->GetTofM();
+ v9 = mu_pl->GetMomentum()->P();
+
+ double y1 = reader->EvaluateMVA( "BDTG-1"); // ML evaluation
+ if (y1 < .72) continue; // Using a cut on the ML output
+ }
+
+ TLorentzVector* P_pl = mu_pl->GetMomentum();
+
+ Double_t mass = mu_pl->GetTofM();
+ Double_t mom = P_pl->P();
+
+ if (type == 0)
+ if (mass < (p0min + p1min * mom + p2min * mom * mom) || mass > (p0max + p1max * mom + p2max * mom * mom)) continue;
+ P1 = *mu_pl->GetMomentum();
+ Plus->Fill();
+ }
+ for (int jPart = 0; jPart < NofMinus; jPart++) {
+ CbmAnaMuonCandidate* mu_mn = (CbmAnaMuonCandidate*) MuMinus->At(jPart);
+ if (type == 0) {
+ if (mu_mn->GetNStsHits() < STShits) continue;
+ if (mu_mn->GetNMuchHits() < MUCHhits) continue;
+ if (mu_mn->GetNTrdHits() < TRDhits) continue;
+ if (mu_mn->GetNTofHits() < TOFhits) continue;
+
+ if (mu_mn->GetChiToVertex() > Vchi2) continue;
+ if (mu_mn->GetChiSts() > STSchi2) continue;
+ if (mu_mn->GetChiMuch() > MUCHchi2) continue;
+
+ v1 = mu_mn->GetNMuchHits();
+ v2 = mu_mn->GetNStsHits();
+ v3 = mu_mn->GetNTrdHits();
+ v4 = mu_mn->GetChiMuch();
+ v5 = mu_mn->GetChiSts();
+ v6 = mu_mn->GetChiToVertex();
+ v7 = mu_mn->GetChiTrd();
+ v8 = mu_mn->GetTofM();
+ v9 = mu_mn->GetMomentum()->P();
+
+ double y1 = reader->EvaluateMVA( "BDTG-1"); // ML evaluation
+ if (y1 < .72) continue; // Using a cut on the ML output
+ }
+ TLorentzVector* P_mn = mu_mn->GetMomentum();
+
+ Double_t mass = mu_mn->GetTofM();
+ Double_t mom = P_mn->P();
+
+ if (type == 0)
+ if (mass < (p0min + p1min * mom + p2min * mom * mom) || mass > (p0max + p1max * mom + p2max * mom * mom)) continue;
+ P2 = *mu_mn->GetMomentum();
+ Minus->Fill();
+ }
+ }
+ f->Close();
+ if (Plus->GetEntries() > 10000) break; // speed up of calculation
+ }
+
+ Plus->SetBranchAddress("P0", &P1(0));
+ Minus->SetBranchAddress("P0", &P2(0));
+
+ Double_t entriesPlus = Plus->GetEntries();
+ Double_t entriesMinus = Minus->GetEntries();
+
+ for (int jPart = 0; jPart < entriesPlus; jPart++) {
+ if (jPart % 1000 == 0) cout << entriesPlus << " ----- " << jPart << endl;
+ Plus->GetEntry(jPart);
+ for (int iPart = 0; iPart < entriesMinus; iPart++) {
+ Minus->GetEntry(iPart);
+ M = P1 + P2;
+ invM_bg->Fill(M.M(), 1. / 10.);
+ }
+ }
+
+ invM_sgn->Scale(1. / NEvents_sgn);
+ invM_full->Add(invM_sgn);
+
+ invM_bg->Scale(1. / NEvents_bg / NEvents_bg);
+
+
+ invM_full->Add(invM_bg);
+
+ Double_t initPar[] = {0.782, 0.02}; // omega mass and sigma for fit initialization
+
+ TF1* fit = new TF1("fullFit", "gaus(0)+pol2(3)", 0.60, .90);
+
+ fit->SetParNames("Const", "Mass", "Sigma", "a0", "a1", "a2");
+ fit->SetNpx(10000);
+ fit->SetLineColor(2);
+ fit->SetParameter(1, (0.60 + .90) / 2);
+ fit->SetParameter(2, (.90 - 0.60) / 4);
+ fit->SetParLimits(2, 0, 1);
+
+ fit->FixParameter(1, initPar[0]);
+ fit->FixParameter(2, initPar[1]);
+ invM_full->Fit("fullFit", "QRN", "", 0.60, .90);
+
+ fit->ReleaseParameter(1);
+ fit->ReleaseParameter(2);
+ invM_full->Fit("fullFit", "QRN", "", 0.60, .90);
+
+ float start = fit->GetParameter(1) - 2 * fit->GetParameter(2);
+ float end = fit->GetParameter(1) + 2 * fit->GetParameter(2);
+
+ cout<<" start "<<start<<" end "<<end<<endl;
+
+ invM_full->Fit("fullFit", "", "", 0.60, .90);
+
+ TF1* fitSgn = new TF1("sgnFit", "gaus", 0.60, .90);
+ for (int iPar = 0; iPar < 3; iPar++)
+ fitSgn->SetParameter(iPar, fit->GetParameter(iPar));
+
+ TF1* fitBg = new TF1("bgFit", "pol2", 0.60, .90);
+ for (int iPar = 3; iPar < 6; iPar++)
+ fitBg->SetParameter(iPar - 3, fit->GetParameter(iPar));
+
+ Double_t min_histo = (invM_sgn->GetXaxis())->GetXmin();
+ Double_t max_histo = (invM_sgn->GetXaxis())->GetXmax();
+ Double_t bin_histo = (max_histo - min_histo) / invM_sgn->GetNbinsX();
+
+ float S = fitSgn->Integral(start, end) / bin_histo;
+ float B = fitBg->Integral(start, end) / bin_histo;
+ float SBratio = S / B;
+ cout<<" signal "<<S<<" Bg "<<B<<" SbyB "<<SBratio<<endl;
+ Double_t eff = S / sgn_mult * 10. * 100.;
+ cout<<"Eff "<<eff<<" % "<<endl;
+ char tmpString[5];
+ sprintf(tmpString, "%4.2f", (S / B));
+
+ TString LegendEntry[6];
+ LegendEntry[0] = "#omega(782)";
+ LegendEntry[1] = "S/B = ";
+ LegendEntry[1] += tmpString;
+ LegendEntry[2] = "#varepsilon = ";
+ sprintf(tmpString, "%2.2f ", eff);
+ LegendEntry[2] += tmpString;
+ LegendEntry[2] += "%";
+ LegendEntry[3] = "#frac{S}{#sqrt{S+B}} = ";
+ sprintf(tmpString, "%1.2e", S / sqrt(S + B));
+ LegendEntry[3] += tmpString;
+ LegendEntry[4] = " ";
+
+ TCanvas* c1;
+ c1 = new TCanvas("c1", "Particles", 0, 0, 1200, 600);
+ c1->UseCurrentStyle();
+ gPad->SetFixedAspectRatio(kFALSE);
+
+ TLegend* legend = legend = new TLegend(0.60, 0.55, .99, .99);
+ legend->SetFillColor(0);
+ legend->SetTextSize(0.05);
+ legend->SetMargin(0.1);
+
+ TGaxis::SetMaxDigits(3);
+
+ invM_full->GetXaxis()->SetTitleSize(0.05);
+ invM_full->GetXaxis()->SetLabelSize(0.04);
+ invM_full->GetXaxis()->SetTitleOffset(0.9);
+ invM_full->GetXaxis()->SetRangeUser(0.6, 1.2);
+
+ invM_full->GetYaxis()->SetTitleSize(0.04);
+ invM_full->GetYaxis()->SetLabelSize(0.04);
+ invM_full->GetYaxis()->SetTitleOffset(1.3);
+
+ invM_full->SetLineColor(kBlack);
+ invM_full->SetMarkerStyle(20);
+ invM_full->SetMarkerColor(kBlack);
+ invM_full->SetMarkerSize(0.7);
+ legend->AddEntry(invM_full, LegendEntry[0].Data(), "l");
+ legend->AddEntry((TObject*) 0, LegendEntry[1].Data(), "");
+ legend->AddEntry((TObject*) 0, LegendEntry[2].Data(), "");
+ legend->AddEntry((TObject*) 0, LegendEntry[3].Data(), "");
+ legend->AddEntry((TObject*) 0, LegendEntry[4].Data(), "");
+
+ if (type == 0) {
+ c1->Divide(2, 1);
+ c1->cd(1);
+
+ invM_full->Draw("pe");
+ legend->Draw();
+
+ c1->cd(2);
+
+ TLatex Tl;
+ Tl.SetTextAlign(12);
+ Tl.SetTextSize(0.04);
+
+ name = "Cuts:";
+ Tl.DrawLatex(0.3, 0.8, name);
+ name.Form("N of STS hits #geq %d", STShits);
+ Tl.DrawLatex(0.3, 0.7, name);
+ name.Form("N of MUCH hits #geq %d", MUCHhits);
+ Tl.DrawLatex(0.3, 0.6, name);
+ name.Form("N of TRD hits #geq %d", TRDhits);
+ Tl.DrawLatex(0.3, 0.5, name);
+ name.Form("N of TOF hits #geq %d", TOFhits);
+ Tl.DrawLatex(0.3, 0.4, name);
+ name.Form("#chi^{2}_{vertex} #leq %1.1f", Vchi2);
+ Tl.DrawLatex(0.3, 0.3, name);
+ name.Form("#chi^{2}_{STS} #leq %1.1f", STSchi2);
+ Tl.DrawLatex(0.3, 0.2, name);
+ name.Form("#chi^{2}_{MUCH} #leq %1.1f", MUCHchi2);
+ Tl.DrawLatex(0.3, 0.1, name);
+ name.Form("%d#sigma cut in TOF", sigmaTOFCut);
+ Tl.DrawLatex(0.3, 0.03, name);
+ }
+ else {
+ invM_full->Draw("pe");
+ legend->Draw();
+ }
+
+ if (type == 0) name.Form("invM_bg_omega_%dgev.cuts_ML_", energy);
+ else
+ name.Form("invM_bg_omega_%dgev.geo_", energy);
+ name += version + ".root";
+
+ TFile* FFF = new TFile(name, "recreate");
+ c1->Write();
+ invM_sgn->Write();
+ invM_bg->Write();
+ invM_full->Write();
+
+ FFF->Close();
+
+ name += ".ps";
+ c1->Print(name);
+}