From 79971c0929f7919f6d959e7f278cfb46d8e757d5 Mon Sep 17 00:00:00 2001 From: Marco van Leeuwen Date: Thu, 22 Jul 2021 14:18:18 +0200 Subject: [PATCH 01/10] Add parametrised response for preshower detector --- src/CMakeLists.txt | 2 ++ src/PreShower.cc | 69 ++++++++++++++++++++++++++++++++++++++++++++++ src/PreShower.hh | 49 ++++++++++++++++++++++++++++++++ 3 files changed, 120 insertions(+) create mode 100644 src/PreShower.cc create mode 100644 src/PreShower.hh diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8a2959b..eef564c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,6 +7,7 @@ set(SOURCES TOFLayer.cc RICHdetector.cc MIDdetector.cc + PreShower.cc ) set(HEADERS @@ -16,6 +17,7 @@ set(HEADERS TOFLayer.hh RICHdetector.hh MIDdetector.hh + PreShower.hh ) get_target_property(DELPHES_INCLUDE_DIRECTORIES diff --git a/src/PreShower.cc b/src/PreShower.cc new file mode 100644 index 0000000..aa161f3 --- /dev/null +++ b/src/PreShower.cc @@ -0,0 +1,69 @@ +/// @author: Roberto Preghenella +/// @email: preghenella@bo.infn.it + +/// @author: Antonio Uras +/// @email: antonio.uras@cern.ch + +/// @author: Marco van Leeuwen +/// @email: marco.van.leeuwen@cern.ch + +#include "PreShower.hh" +#include "TDatabasePDG.h" +#include "THnSparse.h" +#include "TRandom.h" +#include "TDatime.h" +#include "TFile.h" +#include "TVector3.h" +#include "TMath.h" +#include "TAxis.h" + +namespace o2 { + namespace delphes { + + //========================================================================================================== + + bool PreShower::setup() { + + TDatime t; + gRandom->SetSeed(t.GetDate()+t.GetYear()*t.GetHour()*t.GetMinute()*t.GetSecond()); // NB: gRandom is a global pointer ? + for (Int_t iPart = 0; iPart < kNPart; iPart++) { + mMomMin[iPart] = 0.1; + mMomMax[iPart] = 20; + } + + return kTRUE; + + } + + //========================================================================================================== + + bool PreShower::hasPreShower(const Track &track) { + + auto pdg = std::abs(track.PID); + auto part = pidmap[pdg]; + return ((TMath::Abs(track.Eta) < mEtaMax) && (track.P > mMomMin[part])); + + } + + //========================================================================================================== + + bool PreShower::isElectron(const Track &track, int multiplicity=1) { + + auto pdg = std::abs(track.PID); + auto part = pidmap[pdg]; + if (part == kElectron) { + // Parametrisation of preshower detector studies without charge sharing + float eff = 0.8*(1.-exp(1.6*(track.P-0.05))); + return (gRandom->Uniform() < eff); + } + else { + const Float_t misTagProb = 0.001; + return (gRandom->Uniform() < misTagProb); + } + } + + //========================================================================================================== + + } /** namespace delphes **/ + +} /** namespace o2 **/ diff --git a/src/PreShower.hh b/src/PreShower.hh new file mode 100644 index 0000000..0c71de5 --- /dev/null +++ b/src/PreShower.hh @@ -0,0 +1,49 @@ +/// @author: Roberto Preghenella +/// @email: preghenella@bo.infn.it + +/// @author: Antonio Uras +/// @email: antonio.uras@cern.ch + +/// @author: Marco van Leeuwen +/// @email: marco.van.leeuwen@cern.ch + +#ifndef _DelphesO2_PreShower_h_ +#define _DelphesO2_PreShower_h_ + +#include "classes/DelphesClasses.h" +#include "THnSparse.h" +#include "TFile.h" + +#include +using namespace std; + +namespace o2 { + namespace delphes { + + class PreShower { + + public: + PreShower() = default; + ~PreShower() = default; + + enum { kElectron, kMuon, kPion, kKaon, kProton, kNPart }; // primary particles with a non-zero muon PID probability + + bool setup(); + bool hasPreShower(const Track &track); + bool isElectron(const Track &track, int multiplicity); + + protected: + + const double mEtaMax = 1.75; + double mMomMin[kNPart]; + double mMomMax[kNPart]; + const char *partLabel[kNPart] = {"electron","muon","pion","kaon","proton"}; + std::map pidmap = { {11, kElectron}, {13, kMuon}, {211, kPion}, {321, kKaon}, {2212, kProton} }; + + }; + + } /** namespace delphes **/ +} /** namespace o2 **/ + +#endif /** _DelphesO2_MIDLayer_h_ **/ + From d32b13dd1f276a8431463247798dd79610c92ab0 Mon Sep 17 00:00:00 2001 From: Marco van Leeuwen Date: Wed, 28 Jul 2021 22:45:03 +0200 Subject: [PATCH 02/10] Forward tracking parametrisation based on inter/extrapolation; geometry v1.2 --- src/lutWrite.cc | 79 +++++++++++++++++++++++++++++++++++++++------ src/lutWrite.v12.cc | 60 ++++++++++++++++++++++++++++++++++ 2 files changed, 129 insertions(+), 10 deletions(-) create mode 100644 src/lutWrite.v12.cc diff --git a/src/lutWrite.cc b/src/lutWrite.cc index 3972c69..3ce934f 100644 --- a/src/lutWrite.cc +++ b/src/lutWrite.cc @@ -6,6 +6,7 @@ DetectorK fat; void diagonalise(lutEntry_t &lutEntry); +static float etaMaxBarrel = 1.75; bool fatSolve(float *eff, float *covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, int layer = 0, int what = 0, int efftype = 0, int q = 1) @@ -48,8 +49,59 @@ fwdSolve(float *covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000) return true; } +bool +fwdPara(float *covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, float Bfield = 0.5, int layer = 0, int what = 0) +{ + // parametrised forward response; interpolates between FAT at eta = 1.75 and a fixed parametrisation at eta = 4; only diagonal elements + float eff = 0; + float covmbarrel[15] = {0}; + if (fabs(eta) < etaMaxBarrel || fabs(eta) > 4) + return false; + + if (!fatSolve(&eff, covmbarrel, pt, etaMaxBarrel, mass, layer, what)) return false; + + // parametrisation at eta = 4 + double beta = 1./sqrt(1+mass*mass/pt/pt/cosh(eta)/cosh(eta)); + float dca_pos = 2.5e-4/sqrt(3); // 2.5 micron/sqrt(3) + float r0 = 0.5; // layer 0 radius [cm] + float r1 = 1.3; + float r2 = 2.5; + float x0layer = 0.001; // material budget (rad length) per layer + double sigma_alpha = 0.0136/beta/pt*sqrt(x0layer*cosh(eta))*(1+0.038*log(x0layer*cosh(eta))); + double dcaxy_ms = sigma_alpha*r0*sqrt(1+r1*r1/(r2-r0)/(r2-r0)); + double dcaxy2 = dca_pos*dca_pos+dcaxy_ms*dcaxy_ms; + + double dcaz_ms = sigma_alpha*r0*cosh(eta); + double dcaz2 = dca_pos*dca_pos+dcaz_ms*dcaz_ms; + + float Leta = 2.8/sinh(eta)-0.01*r0; // m + double relmomres_pos = 10e-6*pt/0.3/Bfield/Leta/Leta*sqrt(720./15.); + + float relmomres_barrel = sqrt(covmbarrel[14])*pt; + float Router = 1; // m + float relmomres_pos_barrel = 10e-6*pt/0.3/Bfield/Router/Router/sqrt(720./15.); + float relmomres_MS_barrel = sqrt(relmomres_barrel*relmomres_barrel-relmomres_pos_barrel*relmomres_pos_barrel); + + // interpolate MS contrib (rel resolution 0.4 at eta = 4) + float relmomres_MS = 1./beta*0.4*pow(0.4/(beta*relmomres_MS_barrel),(fabs(eta)-4.)/(4.-etaMaxBarrel)); + float momres_tot = pt*sqrt(relmomres_pos*relmomres_pos + relmomres_MS*relmomres_MS); // total absolute mom reso + + // Fill cov matrix diag + for (int i = 0; i < 15; ++i) + covm[i] = 0; + + covm[0] = covmbarrel[0]; + if (dcaxy2 > covm[0]) covm[0] = dcaxy2; + covm[2] = covmbarrel[2]; + if (dcaz2 > covm[2]) covm[2] = dcaz2; + covm[5] = covmbarrel[5]; // sigma^2 sin(phi) + covm[9] = covmbarrel[9]; // sigma^2 tanl + covm[14] = momres_tot*momres_tot/pt/pt/pt/pt; // sigma^2 1/pt + return true; +} + void -lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, int layer = 0, int what = 0, int efftype = 0) +lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, int layer = 0, int what = 0, int efftype = 0, int usePara = 0) { // output file @@ -78,9 +130,9 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, lutHeader.radmap.max = 100.; // eta lutHeader.etamap.log = false; - lutHeader.etamap.nbins = 40; - lutHeader.etamap.min = -2.; - lutHeader.etamap.max = 2.; + lutHeader.etamap.nbins = 80; + lutHeader.etamap.min = -4.; + lutHeader.etamap.max = 4.; // pt lutHeader.ptmap.log = true; lutHeader.ptmap.nbins = 200; @@ -108,7 +160,7 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, for (int ipt = 0; ipt < npt; ++ipt) { lutEntry.pt = lutHeader.ptmap.eval(ipt); lutEntry.valid = true; - if (fabs(eta) < 2.) { + if (fabs(eta) <= etaMaxBarrel) { // full lever arm ends at etaMaxBarrel // printf(" --- fatSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); if (!fatSolve(&lutEntry.eff, lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass, layer, what, efftype, q)) { // printf(" --- fatSolve: error \n"); @@ -120,11 +172,18 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, else { // printf(" --- fwdSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); lutEntry.eff = 1.; - if (!fwdSolve(lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass)) { - // printf(" --- fwdSolve: error \n"); - lutEntry.valid = false; - for (int i = 0; i < 15; ++i) - lutEntry.covm[i] = 0.; + bool retval = true; + if (usePara) { + retval = fwdPara(lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass, field, layer, what); + } + else { + retval = fwdSolve(lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass); + } + if ( !retval) { + // printf(" --- fwdSolve: error \n"); + lutEntry.valid = false; + for (int i = 0; i < 15; ++i) + lutEntry.covm[i] = 0.; } } diagonalise(lutEntry); diff --git a/src/lutWrite.v12.cc b/src/lutWrite.v12.cc new file mode 100644 index 0000000..4e6d27a --- /dev/null +++ b/src/lutWrite.v12.cc @@ -0,0 +1,60 @@ +/// @author: Roberto Preghenella +/// @email: preghenella@bo.infn.it + +#include "lutWrite.cc" + +float scale = 1.; + +void +fatInit_v12(float field = 0.5, float rmin = 100.) +{ + fat.SetBField(field); + // new ideal Pixel properties? + Double_t x0IB = 0.001; + Double_t x0OB = 0.01; + Double_t xrhoIB = 2.3292e-02; // 100 mum Si + Double_t xrhoOB = 2.3292e-01; // 1000 mum Si + + Double_t resRPhiIB = 0.00025; + Double_t resZIB = 0.00025; + Double_t resRPhiOB = 0.00100; + Double_t resZOB = 0.00100; + Double_t eff = 0.98; + fat.AddLayer((char*)"vertex", 0.0, 0, 0); // dummy vertex for matrix calculation + fat.AddLayer((char*)"bpipe0", 0.48 * scale, 0.00042, 2.772e-02); // 150 mum Be + // + fat.AddLayer((char*)"B00", 0.50 * scale, x0IB, xrhoIB, resRPhiIB, resZIB, eff); + fat.AddLayer((char*)"B01", 1.20 * scale, x0IB, xrhoIB, resRPhiIB, resZIB, eff); + fat.AddLayer((char*)"B02", 2.50 * scale, x0IB, xrhoIB, resRPhiIB, resZIB, eff); + // + fat.AddLayer((char*)"bpipe1", 3.7 * scale, 0.0014, 9.24e-02); // 500 mum Be + // + fat.AddLayer((char*)"B03", 3.75 * scale, x0OB, xrhoOB, resRPhiOB, resZOB, eff); + fat.AddLayer((char*)"B04", 7.00 * scale, x0OB, xrhoOB, resRPhiOB, resZOB, eff); + fat.AddLayer((char*)"B05", 12.0 * scale, x0OB, xrhoOB, resRPhiOB, resZOB, eff); + fat.AddLayer((char*)"B06", 20.0 * scale, x0OB, xrhoOB, resRPhiOB, resZOB, eff); + fat.AddLayer((char*)"B07", 30.0 * scale, x0OB, xrhoOB, resRPhiOB, resZOB, eff); + fat.AddLayer((char*)"B08", 45.0 * scale, x0OB, xrhoOB, resRPhiOB, resZOB, eff); + fat.AddLayer((char*)"B09", 60.0 * scale, x0OB, xrhoOB, resRPhiOB, resZOB, eff); + fat.AddLayer((char*)"B10", 80.0 * scale, x0OB, xrhoOB, resRPhiOB, resZOB, eff); + fat.AddLayer((char*)"B11", 100. * scale, x0OB, xrhoOB, resRPhiOB, resZOB, eff); + fat.SetAtLeastHits(4); + fat.SetAtLeastCorr(4); + fat.SetAtLeastFake(0); + // + fat.SetMinRadTrack(rmin); + // + fat.PrintLayout(); +} + +void +lutWrite_v12(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.5, float rmin = 100.) +{ + + // init FAT + fatInit_v12(field, rmin); + // write + lutWrite(filename, pdg, field, 0, 0, 0, 1); // last flag: use parametrisation for forward resolution + +} + From 1248dae7fc500bc5b6616d75e9d8812b2e93c73c Mon Sep 17 00:00:00 2001 From: Marco van Leeuwen Date: Wed, 28 Jul 2021 22:47:45 +0200 Subject: [PATCH 03/10] Plotting macros for forward rapidity and eta dependence --- examples/smearing/draw_etadep.C | 111 +++++++++++++++++++++ examples/smearing/draw_mid_fwd.C | 120 +++++++++++++++++++++++ examples/smearing/lutRead_allres_vseta.C | 106 ++++++++++++++++++++ examples/smearing/lutRead_allres_vspt.C | 106 ++++++++++++++++++++ examples/smearing/lutRead_pt.C | 25 +++-- 5 files changed, 459 insertions(+), 9 deletions(-) create mode 100644 examples/smearing/draw_etadep.C create mode 100644 examples/smearing/draw_mid_fwd.C create mode 100644 examples/smearing/lutRead_allres_vseta.C create mode 100644 examples/smearing/lutRead_allres_vspt.C diff --git a/examples/smearing/draw_etadep.C b/examples/smearing/draw_etadep.C new file mode 100644 index 0000000..4a42777 --- /dev/null +++ b/examples/smearing/draw_etadep.C @@ -0,0 +1,111 @@ +#include "style.C" +#include "lutRead_allres_vseta.C" + +void +draw_etadep() +{ + + style(); + + //std::vector name = {"el", "pi", "ka", "pr"}; + // std::vector title = {"electron", "pion", "kaon", "proton"}; + //std::vector title = {"e", "#pi", "K", "p"}; + + float ptvals[] = {1, 5, 10, 20}; + Color_t colors[] = {kRed+1, kRed+1, kRed+1, kRed+1}; + Color_t linestyles[] = {1, 2, 3, 4}; + + int nPtVals = 4; + + TLatex latex; + latex.SetTextAlign(33); + + auto c1 = new TCanvas("c1", "c1: pt reso", 800, 600); + c1->SetLogy(); + c1->DrawFrame(0, 1., 4, 200., ";#eta;momentum resolution (%)"); + + TLegend *leg = new TLegend(0.18,0.6,0.3,0.85); + leg->SetBorderSize(0); + auto c2 = new TCanvas("c2", "c2: dca xy", 800, 600); + c2->SetLogy(); + c2->DrawFrame(0, 1., 4, 200., ";#eta;dca_{xy} resolution (#mum)"); + + auto c3 = new TCanvas("c3", "c3: dca z", 800, 600); + c3->SetLogy(); + c3->DrawFrame(0, 1., 4, 200., ";#eta;dca_{z} resolution (#mum)"); + + auto c4 = new TCanvas("c4", "c4: sin phi", 800, 600); + c4->SetLogy(); + c4->DrawFrame(0, 1e-5, 4, 1e-2, ";#eta;sin(#phi) resolution"); + + auto c5 = new TCanvas("c5", "c5: tan lambda", 800, 600); + c5->SetLogy(); + c5->DrawFrame(0, 1e-5, 4, 1e-2, ";#eta;tan(#lambda) resolution"); + + + for (int i = 0; i < nPtVals; ++i) { + + /* + auto c = new TCanvas((std::string("c") + name[i]).c_str(), + (std::string("c") + name[i]).c_str(), + 800, 800); + c->SetLogx(); + c->SetLogy(); + c->DrawFrame(1.e-2, 1., 100., 100., ";#it{p}_{T} (GeV/#it{c});momentum resolution (%)"); + latex.DrawLatexNDC(0.9, 0.9, title[i].c_str()); + */ + + //auto g2a = lutRead_pt((std::string("lutCovm.") + name[i] + std::string(".2kG.20cm.default.dat")).c_str()); + auto ptres = new TGraph(); + auto dca_xy = new TGraph(); + auto dca_z = new TGraph(); + auto sinp = new TGraph(); + auto tanl = new TGraph(); + lutRead_allres_vseta("luts/lutCovm.v12.dat",dca_xy,dca_z,sinp,tanl,ptres,ptvals[i]); + ptres->SetLineColor(colors[i]); + ptres->SetLineStyle(i+1); + ptres->SetLineWidth(3); + leg->AddEntry(ptres,Form("p_{T} = %.1f GeV",ptvals[i])); + c1->cd(); + ptres->Draw("samel"); + + dca_xy->SetLineColor(colors[i]); + dca_xy->SetLineStyle(i+1); + dca_xy->SetLineWidth(3); + c2->cd(); + dca_xy->Draw("samel"); + + dca_z->SetLineColor(colors[i]); + dca_z->SetLineStyle(i+1); + dca_z->SetLineWidth(3); + c3->cd(); + dca_z->Draw("samel"); + + sinp->SetLineColor(colors[i]); + sinp->SetLineStyle(i+1); + sinp->SetLineWidth(3); + c4->cd(); + sinp->Draw("samel"); + + tanl->SetLineColor(colors[i]); + tanl->SetLineStyle(i+1); + tanl->SetLineWidth(3); + c5->cd(); + tanl->Draw("samel"); + } + c1->cd(); + leg->Draw(); + c1->Print("ptres_etadep.pdf"); + c2->cd(); + leg->Draw(); + c2->Print("dcaxyres_etadep.pdf"); + c3->cd(); + leg->Draw(); + c3->Print("dcazres_etadep.pdf"); + c4->cd(); + leg->Draw(); + c4->Print("sinphires_etadep.pdf"); + c5->cd(); + leg->Draw(); + c5->Print("tanlres_etadep.pdf"); +} diff --git a/examples/smearing/draw_mid_fwd.C b/examples/smearing/draw_mid_fwd.C new file mode 100644 index 0000000..e2becf0 --- /dev/null +++ b/examples/smearing/draw_mid_fwd.C @@ -0,0 +1,120 @@ +#include "style.C" +#include "lutRead_allres_vspt.C" + +void +draw_mid_fwd() +{ + + style(); + + //std::vector name = {"el", "pi", "ka", "pr"}; + // std::vector title = {"electron", "pion", "kaon", "proton"}; + //std::vector title = {"e", "#pi", "K", "p"}; + + float etavals[] = {0, 1, 2, 3, 4}; + Color_t colors[] = {kRed+1, kRed+1, kBlue, kBlue, kBlue}; + Color_t linestyles[] = {1, 2, 3, 4, 5}; + + int nEtaVals = 5; + + TLatex latex; + latex.SetTextAlign(33); + + auto c1 = new TCanvas("c1", "c1: pt reso", 800, 600); + //cc->Divide(2, 2); + c1->SetLogx(); + c1->SetLogy(); + c1->DrawFrame(1.e-2, 1., 100., 200., ";#it{p}_{T} (GeV/#it{c});momentum resolution (%)"); + + TLegend *leg = new TLegend(0.18,0.18,0.3,0.35); + leg->SetBorderSize(0); + auto c2 = new TCanvas("c2", "c2: dca xy", 800, 600); + c2->SetLogx(); + c2->SetLogy(); + c2->DrawFrame(1.e-2, 1., 100., 200., ";#it{p}_{T} (GeV/#it{c});dca_{xy} resolution (#mum)"); + + auto c3 = new TCanvas("c3", "c3: dca z", 800, 600); + c3->SetLogx(); + c3->SetLogy(); + c3->DrawFrame(1.e-2, 1., 100., 200., ";#it{p}_{T} (GeV/#it{c});dca_{z} resolution (#mum)"); + + auto c4 = new TCanvas("c4", "c4: sin phi", 800, 600); + c4->SetLogx(); + c4->SetLogy(); + c4->DrawFrame(1.e-2, 0.00001, 100., 5e-2, ";#it{p}_{T} (GeV/#it{c});sin(#phi) resolution"); + + auto c5 = new TCanvas("c5", "c5: tan lambda", 800, 600); + c5->SetLogx(); + c5->SetLogy(); + c5->DrawFrame(1.e-2, 0.00001, 100., 5e-2, ";#it{p}_{T} (GeV/#it{c});tan(#lambda) resolution"); + + + for (int i = 0; i < nEtaVals; ++i) { + + /* + auto c = new TCanvas((std::string("c") + name[i]).c_str(), + (std::string("c") + name[i]).c_str(), + 800, 800); + c->SetLogx(); + c->SetLogy(); + c->DrawFrame(1.e-2, 1., 100., 100., ";#it{p}_{T} (GeV/#it{c});momentum resolution (%)"); + latex.DrawLatexNDC(0.9, 0.9, title[i].c_str()); + */ + + //auto g2a = lutRead_pt((std::string("lutCovm.") + name[i] + std::string(".2kG.20cm.default.dat")).c_str()); + auto ptres = new TGraph(); + auto dca_xy = new TGraph(); + auto dca_z = new TGraph(); + auto sinp = new TGraph(); + auto tanl = new TGraph(); + lutRead_allres_vspt("luts/lutCovm.v12.dat",dca_xy,dca_z,sinp,tanl,ptres,etavals[i]); + ptres->SetLineColor(colors[i]); + ptres->SetLineStyle(i+1); + ptres->SetLineWidth(3); + leg->AddEntry(ptres,Form("#eta = %.1f",etavals[i])); + c1->cd(); + ptres->Draw("samel"); + + dca_xy->SetLineColor(colors[i]); + dca_xy->SetLineStyle(i+1); + dca_xy->SetLineWidth(3); + c2->cd(); + dca_xy->Draw("samel"); + + dca_z->SetLineColor(colors[i]); + dca_z->SetLineStyle(i+1); + dca_z->SetLineWidth(3); + c3->cd(); + dca_z->Draw("samel"); + + sinp->SetLineColor(colors[i]); + sinp->SetLineStyle(i+1); + sinp->SetLineWidth(3); + c4->cd(); + sinp->Draw("samel"); + + tanl->SetLineColor(colors[i]); + tanl->SetLineStyle(i+1); + tanl->SetLineWidth(3); + c5->cd(); + tanl->Draw("samel"); + } + c1->cd(); + leg->Draw(); + c1->Print("ptres_etabins.pdf"); + c2->cd(); + leg->Draw(); + c2->Print("dcaxyres_etabins.pdf"); + c3->cd(); + leg->Draw(); + c3->Print("dcazres_etabins.pdf"); + c4->cd(); + leg->Draw(); + c4->Print("sinphires_etabins.pdf"); + c5->cd(); + leg->Draw(); + c5->Print("tanlres_etabins.pdf"); + + //cc->SaveAs("draw_pt_mid_fwd.png"); + +} diff --git a/examples/smearing/lutRead_allres_vseta.C b/examples/smearing/lutRead_allres_vseta.C new file mode 100644 index 0000000..565c4d9 --- /dev/null +++ b/examples/smearing/lutRead_allres_vseta.C @@ -0,0 +1,106 @@ +#include "lutCovm.hh" + +void lutRead_allres_vseta(const char *filename = "lutCovm.dat", + TGraph *dca_xy = 0, TGraph *dca_z = 0, TGraph *sinp = 0, TGraph *tanl = 0, TGraph* ptres = 0, double pt = 1.) +{ + + // input file + ifstream lutFile(filename, std::ofstream::binary); + + // read header + lutHeader_t lutHeader; + lutFile.read(reinterpret_cast(&lutHeader), sizeof(lutHeader)); + lutHeader.print(); + + // entries + const int nnch = lutHeader.nchmap.nbins; + const int nrad = lutHeader.radmap.nbins; + const int neta = lutHeader.etamap.nbins; + const int npt = lutHeader.ptmap.nbins; + + auto nchbin = lutHeader.nchmap.find(0.); + auto radbin = lutHeader.radmap.find(0.); + auto ptbin = lutHeader.ptmap.find(pt); + + lutEntry_t lutTable[neta]; + + // read entries + for (int inch = 0; inch < nnch; ++inch) { + for (int irad = 0; irad < nrad; ++irad) { + for (int ieta = 0; ieta < neta; ++ieta) { + for (int ipt = 0; ipt < npt; ++ipt) { + if (inch==nchbin && irad==radbin && ipt == ptbin) { + lutFile.read(reinterpret_cast(&lutTable[ieta]), sizeof(lutEntry_t)); + //lutTable[ieta].print(); + } + else { + lutEntry_t dummy; + lutFile.read(reinterpret_cast(&dummy), sizeof(lutEntry_t)); + } + } + } + } + } + + lutFile.close(); + + if (ptres) { + ptres->SetName(filename); + ptres->SetTitle(filename); + ptres->GetXaxis()->SetTitle("#eta"); + ptres->GetYaxis()->SetTitle("momentum resolution (%)"); + } + + if (dca_xy) { + dca_xy->SetName(filename); + dca_xy->SetTitle(filename); + dca_xy->GetXaxis()->SetTitle("#eta"); + dca_xy->GetYaxis()->SetTitle("#sigma(d_{xy}) (#mum)"); + } + + if (dca_z) { + dca_z->SetName(filename); + dca_z->SetTitle(filename); + dca_z->GetXaxis()->SetTitle("#eta"); + dca_z->GetYaxis()->SetTitle("#sigma(d_{z}) (#mum)"); + } + + if (sinp) { + sinp->SetName(filename); + sinp->SetTitle(filename); + sinp->GetXaxis()->SetTitle("#eta"); + sinp->GetYaxis()->SetTitle("#sigma(sin(#phi))"); + } + + if (tanl) { + tanl->SetName(filename); + tanl->SetTitle(filename); + tanl->GetXaxis()->SetTitle("#eta"); + tanl->GetYaxis()->SetTitle("#sigma(tam(#lambda))"); + } + for (int ieta = 0; ieta < neta/2; ++ieta) { + auto lutEntry = &lutTable[neta/2+ieta]; + if (!lutEntry->valid) continue; + auto cen = lutEntry->eta; + if (ptres) { + auto val = sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.; + ptres->SetPoint(ptres->GetN(), cen, val); + } + if (dca_xy) { + auto val = sqrt(lutEntry->covm[0]) * 1e4; + dca_xy->SetPoint(dca_xy->GetN(), cen, val); + } + if (dca_z) { + auto val = sqrt(lutEntry->covm[2]) * 1e4; + dca_z->SetPoint(dca_z->GetN(), cen, val); + } + if (sinp) { + auto val = sqrt(lutEntry->covm[5]); + sinp->SetPoint(sinp->GetN(), cen, val); + } + if (tanl) { + auto val = sqrt(lutEntry->covm[9]); + tanl->SetPoint(tanl->GetN(), cen, val); + } + } +} diff --git a/examples/smearing/lutRead_allres_vspt.C b/examples/smearing/lutRead_allres_vspt.C new file mode 100644 index 0000000..0cb8776 --- /dev/null +++ b/examples/smearing/lutRead_allres_vspt.C @@ -0,0 +1,106 @@ +#include "lutCovm.hh" + +void lutRead_allres_vspt(const char *filename = "lutCovm.dat", + TGraph *dca_xy = 0, TGraph *dca_z = 0, TGraph *sinp = 0, TGraph *tanl = 0, TGraph* ptres = 0, double eta = 0.) +{ + + // input file + ifstream lutFile(filename, std::ofstream::binary); + + // read header + lutHeader_t lutHeader; + lutFile.read(reinterpret_cast(&lutHeader), sizeof(lutHeader)); + lutHeader.print(); + + // entries + const int nnch = lutHeader.nchmap.nbins; + const int nrad = lutHeader.radmap.nbins; + const int neta = lutHeader.etamap.nbins; + const int npt = lutHeader.ptmap.nbins; + + auto nchbin = lutHeader.nchmap.find(0.); + auto radbin = lutHeader.radmap.find(0.); + auto etabin = lutHeader.etamap.find(eta); + + lutEntry_t lutTable[npt]; + + // read entries + for (int inch = 0; inch < nnch; ++inch) { + for (int irad = 0; irad < nrad; ++irad) { + for (int ieta = 0; ieta < neta; ++ieta) { + for (int ipt = 0; ipt < npt; ++ipt) { + if (inch==nchbin && irad==radbin && ieta == etabin) { + lutFile.read(reinterpret_cast(&lutTable[ipt]), sizeof(lutEntry_t)); + //lutTable[ipt].print(); + } + else { + lutEntry_t dummy; + lutFile.read(reinterpret_cast(&dummy), sizeof(lutEntry_t)); + } + } + } + } + } + + lutFile.close(); + + if (ptres) { + ptres->SetName(filename); + ptres->SetTitle(filename); + ptres->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})"); + ptres->GetYaxis()->SetTitle("momentum resolution (%)"); + } + + if (dca_xy) { + dca_xy->SetName(filename); + dca_xy->SetTitle(filename); + dca_xy->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})"); + dca_xy->GetYaxis()->SetTitle("#sigma(d_{xy}) (#mum)"); + } + + if (dca_z) { + dca_z->SetName(filename); + dca_z->SetTitle(filename); + dca_z->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})"); + dca_z->GetYaxis()->SetTitle("#sigma(d_{z}) (#mum)"); + } + + if (sinp) { + sinp->SetName(filename); + sinp->SetTitle(filename); + sinp->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})"); + sinp->GetYaxis()->SetTitle("#sigma(sin(#phi))"); + } + + if (tanl) { + tanl->SetName(filename); + tanl->SetTitle(filename); + tanl->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})"); + tanl->GetYaxis()->SetTitle("#sigma(tam(#lambda))"); + } + for (int ipt = 0; ipt < npt; ++ipt) { + auto lutEntry = &lutTable[ipt]; + if (!lutEntry->valid) continue; + auto cen = lutEntry->pt; + if (ptres) { + auto val = sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.; + ptres->SetPoint(ptres->GetN(), cen, val); + } + if (dca_xy) { + auto val = sqrt(lutEntry->covm[0]) * 1e4; + dca_xy->SetPoint(dca_xy->GetN(), cen, val); + } + if (dca_z) { + auto val = sqrt(lutEntry->covm[2]) * 1e4; + dca_z->SetPoint(dca_z->GetN(), cen, val); + } + if (sinp) { + auto val = sqrt(lutEntry->covm[5]); + sinp->SetPoint(sinp->GetN(), cen, val); + } + if (tanl) { + auto val = sqrt(lutEntry->covm[9]); + tanl->SetPoint(tanl->GetN(), cen, val); + } + } +} diff --git a/examples/smearing/lutRead_pt.C b/examples/smearing/lutRead_pt.C index c57acb1..43a871c 100644 --- a/examples/smearing/lutRead_pt.C +++ b/examples/smearing/lutRead_pt.C @@ -11,39 +11,46 @@ lutRead_pt(const char *filename = "lutCovm.dat", double eta = 0.) lutHeader_t lutHeader; lutFile.read(reinterpret_cast(&lutHeader), sizeof(lutHeader)); lutHeader.print(); - + cout << "header done" << endl; // entries const int nnch = lutHeader.nchmap.nbins; const int nrad = lutHeader.radmap.nbins; const int neta = lutHeader.etamap.nbins; const int npt = lutHeader.ptmap.nbins; - lutEntry_t lutTable[nnch][nrad][neta][npt]; + + auto nchbin = lutHeader.nchmap.find(0.); + auto radbin = lutHeader.radmap.find(0.); + auto etabin = lutHeader.etamap.find(eta); + + lutEntry_t lutTable[npt]; // read entries for (int inch = 0; inch < nnch; ++inch) { for (int irad = 0; irad < nrad; ++irad) { for (int ieta = 0; ieta < neta; ++ieta) { for (int ipt = 0; ipt < npt; ++ipt) { - lutFile.read(reinterpret_cast(&lutTable[inch][irad][ieta][ipt]), sizeof(lutEntry_t)); - // lutTable[inch][irad][ieta][ipt].print(); + if (inch==nchbin && irad==radbin && ieta == etabin) { + lutFile.read(reinterpret_cast(&lutTable[ipt]), sizeof(lutEntry_t)); + //lutTable[ipt].print(); + } + else { + lutEntry_t dummy; + lutFile.read(reinterpret_cast(&dummy), sizeof(lutEntry_t)); + } } } } } lutFile.close(); - // create graph of pt resolution at eta = 0 - auto inch = lutHeader.nchmap.find(0.); - auto irad = lutHeader.radmap.find(0.); - auto ieta = lutHeader.etamap.find(eta); auto gpt = new TGraph(); gpt->SetName(filename); gpt->SetTitle(filename); gpt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})"); gpt->GetYaxis()->SetTitle("momentum resolution (%)"); for (int ipt = 0; ipt < npt; ++ipt) { - auto lutEntry = &lutTable[inch][irad][ieta][ipt]; + auto lutEntry = &lutTable[ipt]; if (!lutEntry->valid) continue; auto cen = lutEntry->pt; auto val = sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.; From 81358fb033f3ac3e6a75dc201922b1c84475e1b2 Mon Sep 17 00:00:00 2001 From: Marco van Leeuwen Date: Fri, 30 Jul 2021 21:32:48 +0200 Subject: [PATCH 04/10] Include Bfield dependence for MS term at large eta --- src/lutWrite.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/lutWrite.cc b/src/lutWrite.cc index 3ce934f..f83f281 100644 --- a/src/lutWrite.cc +++ b/src/lutWrite.cc @@ -83,7 +83,8 @@ fwdPara(float *covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, f float relmomres_MS_barrel = sqrt(relmomres_barrel*relmomres_barrel-relmomres_pos_barrel*relmomres_pos_barrel); // interpolate MS contrib (rel resolution 0.4 at eta = 4) - float relmomres_MS = 1./beta*0.4*pow(0.4/(beta*relmomres_MS_barrel),(fabs(eta)-4.)/(4.-etaMaxBarrel)); + float relmomres_MS_eta4 = 0.4/beta*0.5/Bfield; + float relmomres_MS = relmomres_MS_eta4*pow(relmomres_MS_eta4/relmomres_MS_barrel,(fabs(eta)-4.)/(4.-etaMaxBarrel)); float momres_tot = pt*sqrt(relmomres_pos*relmomres_pos + relmomres_MS*relmomres_MS); // total absolute mom reso // Fill cov matrix diag From d66fba6c76708e504d7fd7793281d106c8c3e54e Mon Sep 17 00:00:00 2001 From: Roberto Preghenella Date: Tue, 27 Jul 2021 17:24:09 +0200 Subject: [PATCH 05/10] simple number formatting update --- src/DetectorK/DetectorK.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/DetectorK/DetectorK.cxx b/src/DetectorK/DetectorK.cxx index a1ce3f9..57aee34 100644 --- a/src/DetectorK/DetectorK.cxx +++ b/src/DetectorK/DetectorK.cxx @@ -461,11 +461,11 @@ void DetectorK::PrintLayout(Bool_t full) { if (tmp->phiRes==RIDICULOUS) printf(" - "); else - printf("%3.0f ",tmp->phiRes*10000); + printf("%3.2f ",tmp->phiRes*10000); if (tmp->zRes==RIDICULOUS) printf(" -"); else - printf("%3.0f",tmp->zRes*10000); + printf("%3.2f",tmp->zRes*10000); if (tmp->zRes==RIDICULOUS) printf("\t -\n"); From 392b2498e2c1aa6255ebcf146f03c4fde4cf3aec Mon Sep 17 00:00:00 2001 From: Roberto Preghenella Date: Sat, 31 Jul 2021 15:58:01 +0200 Subject: [PATCH 06/10] Get ready for better LUT handling --- src/lutCovm.hh | 6 ++++-- src/lutWrite.cc | 54 ++++++++++++++++++------------------------------- 2 files changed, 24 insertions(+), 36 deletions(-) diff --git a/src/lutCovm.hh b/src/lutCovm.hh index 91af202..8cb8c6a 100644 --- a/src/lutCovm.hh +++ b/src/lutCovm.hh @@ -2,7 +2,7 @@ /// @email: preghenella@bo.infn.it #pragma once -#define LUTCOVM_VERSION 20210611 +#define LUTCOVM_VERSION 20210731 struct map_t { int nbins = 1; @@ -55,7 +55,9 @@ struct lutEntry_t { float eta = 0.; float pt = 0.; bool valid = false; - float eff = 0.; + float eff = 0.; + float itof = 0.; + float otof = 0.; float covm[15] = {0.}; float eigval[5] = {0.}; float eigvec[5][5] = {0.}; diff --git a/src/lutWrite.cc b/src/lutWrite.cc index 47e52da..2fe924f 100644 --- a/src/lutWrite.cc +++ b/src/lutWrite.cc @@ -8,39 +8,23 @@ DetectorK fat; void diagonalise(lutEntry_t &lutEntry); bool -fatSolve(float *eff, float *covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, int layer = 0, int what = 0, int efftype = 0, int q = 1) +fatSolve(lutEntry_t &lutEntry, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, int itof = 0, int otof = 0, int q = 1) { - if (q > 1) { - mass = -mass; - } - TrackSol tr(1, pt, eta, q, mass); - bool retval = fat.SolveTrack(tr); - if (!retval) return false; + lutEntry.valid = false; - switch (efftype) { - case 0: - *eff = fat.GetEfficiency(); - break; - case 1: - *eff = fat.GetLayerEfficiency(layer); - break; - } + // solve track + if (q > 1) mass = -mass; + TrackSol tr(1, pt, eta, q, mass); + if (!fat.SolveTrack(tr)) return false; + AliExternalTrackParam *trPtr = (AliExternalTrackParam*)tr.fTrackCmb.At(0); + if (!trPtr) return false; - AliExternalTrackParam *trPtr = nullptr; - switch (what) { - case 0: - trPtr = (AliExternalTrackParam*)tr.fTrackCmb[layer]; - break; - case 1: - trPtr = (AliExternalTrackParam*)tr.fTrackOutB[layer]; - break; - case 2: - trPtr = (AliExternalTrackParam*)tr.fTrackOutA[layer]; - break; - } + lutEntry.valid = true; + lutEntry.eff = fat.GetGoodHitProb(0); + lutEntry.itof = fat.GetGoodHitProb(itof); + lutEntry.otof = fat.GetGoodHitProb(otof); + for (int i = 0; i < 15; ++i) lutEntry.covm[i] = trPtr->GetCovariance()[i]; - for (int i = 0; i < 15; ++i) - covm[i] = trPtr->GetCovariance()[i]; return true; } @@ -52,7 +36,7 @@ fwdSolve(float *covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000) } void -lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, int layer = 0, int what = 0, int efftype = 0) +lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, int itof = 0, int otof = 0) { // output file @@ -116,10 +100,11 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, lutEntry.pt = lutHeader.ptmap.eval(ipt); lutEntry.valid = true; if (fabs(eta) < 2.) { - // printf(" --- fatSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); - if (!fatSolve(&lutEntry.eff, lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass, layer, what, efftype, q)) { - // printf(" --- fatSolve: error \n"); - lutEntry.valid = false; + // printf(" --- fatSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); + if (!fatSolve(lutEntry, lutEntry.pt, lutEntry.eta, lutHeader.mass, itof, otof, q)) { + // printf(" --- [ERROR] fatSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); + lutEntry.valid = false; + lutEntry.eff = 0.; for (int i = 0; i < 15; ++i) lutEntry.covm[i] = 0.; } @@ -130,6 +115,7 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, if (!fwdSolve(lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass)) { // printf(" --- fwdSolve: error \n"); lutEntry.valid = false; + lutEntry.eff = 0.; for (int i = 0; i < 15; ++i) lutEntry.covm[i] = 0.; } From 57dd6cdd20f766ef061baa939ff225a7d117424a Mon Sep 17 00:00:00 2001 From: Roberto Preghenella Date: Sat, 31 Jul 2021 15:59:02 +0200 Subject: [PATCH 07/10] Universal lutRead macro --- examples/smearing/lutRead.C | 59 +++++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) create mode 100644 examples/smearing/lutRead.C diff --git a/examples/smearing/lutRead.C b/examples/smearing/lutRead.C new file mode 100644 index 0000000..d43fee4 --- /dev/null +++ b/examples/smearing/lutRead.C @@ -0,0 +1,59 @@ +enum EWhat { + kEfficiency, + kEfficiencyInnerTOF, + kEfficiencyOuterTOF, + kPtResolution, + kRPhiResolution, + kZResolution +}; + +enum EVs { + kNch, + kEta, + kPt +}; + +TGraph * +lutRead(int pdg, const char *filename, int what, int vs, float nch = 0., float radius = 0., float eta = 0., float pt = 0.) +{ + o2::delphes::TrackSmearer smearer; + smearer.loadTable(pdg, filename); + auto lutHeader = smearer.getLUTHeader(pdg); + map_t lutMap; + if (vs == kNch) lutMap = lutHeader->nchmap; + if (vs == kEta) lutMap = lutHeader->etamap; + if (vs == kPt) lutMap = lutHeader->ptmap; + auto nbins = lutMap.nbins; + auto g = new TGraph(); + + bool canBeInvalid = true; + for (int i = 0; i < nbins; ++i) { + if (vs == kNch) nch = lutMap.eval(i); + if (vs == kEta) eta = lutMap.eval(i); + if (vs == kPt) pt = lutMap.eval(i); + auto lutEntry = smearer.getLUTEntry(pdg, nch, radius, eta , pt); + if (!lutEntry->valid || lutEntry->eff == 0.) { + if (!canBeInvalid) std::cout << " --- warning: it cannot be invalid \n" << std::endl; + continue; + } + canBeInvalid = false; + + double cen = 0.; + if (vs == kNch) cen = lutEntry->nch; + if (vs == kEta) cen = lutEntry->eta; + if (vs == kPt) cen = lutEntry->pt; + double val = 0.; + if (what == kEfficiency) val = lutEntry->eff * 100.; // efficiency (%) + if (what == kEfficiencyInnerTOF) val = lutEntry->itof * 100.; // efficiency (%) + if (what == kEfficiencyOuterTOF) val = lutEntry->otof * 100.; // efficiency (%) + if (what == kPtResolution) val = sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.; // pt resolution (%) + if (what == kRPhiResolution) val = sqrt(lutEntry->covm[0]) * 1.e4; // rphi resolution (um) + if (what == kZResolution) val = sqrt(lutEntry->covm[1]) * 1.e4; // z resolution (um) + if (val < 0.) continue; + g->SetPoint(g->GetN(), cen, val); + } + + return g; + +} + From 85e1502c4c4862462c9e8374dc0154a0c49d10a8 Mon Sep 17 00:00:00 2001 From: Roberto Preghenella Date: Sat, 31 Jul 2021 15:52:52 +0200 Subject: [PATCH 08/10] Get ready for better LUT handling --- src/DetectorK/DetectorK.cxx | 262 +++++++++++----------- src/DetectorK/DetectorK.h | 13 +- src/DetectorK/HistoManager.cxx | 383 +++++++++++++++++++++++++++++++++ src/DetectorK/HistoManager.h | 74 +++++++ 4 files changed, 605 insertions(+), 127 deletions(-) create mode 100644 src/DetectorK/HistoManager.cxx create mode 100644 src/DetectorK/HistoManager.h diff --git a/src/DetectorK/DetectorK.cxx b/src/DetectorK/DetectorK.cxx index 57aee34..54d42b5 100644 --- a/src/DetectorK/DetectorK.cxx +++ b/src/DetectorK/DetectorK.cxx @@ -27,6 +27,7 @@ Bool_t DetectorK::verboseR=0; #define RIDICULOUS 999999 // A ridiculously large resolution (cm) to flag a dead detector +#define xrhosteps 100 // steps for dEdx correction #define Luminosity 1.e27 // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 ) #define SigmaD 6.0 // Size of the interaction diamond (cm) (LHC = 6.0 cm) #define dNdEtaMinB 1//950//660//950 // Multiplicity per unit Eta (AuAu MinBias = 170, Central = 700) @@ -103,7 +104,7 @@ DetectorK::DetectorK() SetMaxSnp(); } -DetectorK::DetectorK(const char *name, const char *title) +DetectorK::DetectorK(char *name, char *title) : TNamed(name,title), fNumberOfLayers(0), fNumberOfActiveLayers(0), @@ -461,11 +462,11 @@ void DetectorK::PrintLayout(Bool_t full) { if (tmp->phiRes==RIDICULOUS) printf(" - "); else - printf("%3.2f ",tmp->phiRes*10000); + printf("%3.0f ",tmp->phiRes*10000); if (tmp->zRes==RIDICULOUS) printf(" -"); else - printf("%3.2f",tmp->zRes*10000); + printf("%3.0f",tmp->zRes*10000); if (tmp->zRes==RIDICULOUS) printf("\t -\n"); @@ -915,8 +916,8 @@ void DetectorK::SolveViaBilloir(Double_t selPt, double ptmin) { if (lr->xrho>0) { // correct in small steps bool elossOK = kTRUE; - for (int ise=10;ise--;) { - if (!probTrLast.CorrectForMeanMaterial(0, -lr->xrho/10, fParticleMass , kTRUE)) {elossOK = kFALSE; break;} + for (int ise=xrhosteps;ise--;) { + if (!probTrLast.CorrectForMeanMaterial(0, -lr->xrho/xrhosteps, fParticleMass , kTRUE)) {elossOK = kFALSE; break;} } if (!elossOK) break; } @@ -929,14 +930,15 @@ void DetectorK::SolveViaBilloir(Double_t selPt, double ptmin) { lastReached = il; prepLrOK[il] = 1.; // flag successfully passed layer } - // if ( ((CylLayerK*)fLayers.At(lastReached))->radius < fMinRadTrack) continue; + // if ( ((CylLayerK*)fLayers.At(lastReached))->radius < fMinRadTrack) continue; if (!PropagateToR(&probTr,probTr.GetX() + kTrackingMargin,bGauss,1)) continue; // if (probTr.GetX()radius + kTrackingMargin,bGauss,1)) continue; //if (!probTr.PropagateTo(last->radius,bGauss)) continue; // reset cov.matrix @@ -1032,8 +1034,8 @@ void DetectorK::SolveViaBilloir(Double_t selPt, double ptmin) { exit(1); } if (layer->xrho>0) { // correct in small steps - for (int ise=10;ise--;) { - if (!probTr.CorrectForMeanMaterial(0, layer->xrho/10, fParticleMass , kTRUE)) { + for (int ise=xrhosteps;ise--;) { + if (!probTr.CorrectForMeanMaterial(0, layer->xrho/xrhosteps, fParticleMass , kTRUE)) { printf("Failed to apply material correction, xrho=%.4f\n",layer->xrho); probTr.Print(); exit(1); @@ -1245,8 +1247,8 @@ void DetectorK::SolveViaBilloir(Double_t selPt, double ptmin) { exit(1); } if (layer->xrho>0) { // correct in small steps - for (int ise=10;ise--;) { - if (!probTr.CorrectForMeanMaterial(0, -layer->xrho/10, fParticleMass , kTRUE)) { + for (int ise=xrhosteps;ise--;) { + if (!probTr.CorrectForMeanMaterial(0, -layer->xrho/xrhosteps, fParticleMass , kTRUE)) { printf("Failed to apply material correction, xrho=%.4f\n",-layer->xrho); probTr.Print(); exit(1); @@ -1332,7 +1334,7 @@ void DetectorK::SolveViaBilloir(Double_t selPt, double ptmin) { } // Set Detector-Efficiency Storage area to unity fEfficiency[i] = 1.0 ; - + // print out and efficiency calculation iLayActive=0; if (print == 1 && fTransMomenta[i] >= selPt && printOnce == 1) printf("\n Combined propagation estimates\n"); @@ -1424,7 +1426,12 @@ Bool_t DetectorK::SolveTrack(TrackSol& ts) { double etaTr = ts.fEta; double mass = ts.fMass; double charge = ts.fCharge; - + + // reset good hit probability + for (int i = 0; i < kMaxNumberOfDetectors; ++i) + fGoodHitProb[i] = -1.; + fGoodHitProb[0] = 1.; // we use layer zero to accumulate + if (ptTr<0) { printf("Input track is not initialized"); return kFALSE; @@ -1440,7 +1447,6 @@ Bool_t DetectorK::SolveTrack(TrackSol& ts) { TClonesArray &saveParOutwardA = ts.fTrackOutA; TClonesArray &saveParComb = ts.fTrackCmb; - // Calculate track parameters using Billoirs method of matrices Double_t pt,lambda; // CylLayerK *last = (CylLayerK*) fLayers.At((fLayers.GetEntries()-1)); @@ -1482,39 +1488,51 @@ Bool_t DetectorK::SolveTrack(TrackSol& ts) { // // find max layer this track can reach double rmx = (TMath::Abs(fBField)>1e-5) ? pt*100./(0.3*TMath::Abs(fBField)) : 9999; - if (2*rmx-5. < minRad && minRad>0) { + // if (2*rmx-5. < minRad && minRad>0) { + if ( minRad/(2.*rmx)>fMaxSnp-0.01 && minRad>0) { // printf("Track of pt=%.3f cannot be tracked to min. r=%f\n",pt,minRad); return kFALSE; } - Int_t lastActiveLayer = -1; + Int_t lastActiveLayer = -1, lastReachedLayer = -1; for (Int_t j=fLayers.GetEntries(); j--;) { CylLayerK *l = (CylLayerK*) fLayers.At(j); - // printf("at lr %d r: %f vs %f, pt:%f\n",j,l->radius, 2*rmx-2.*kTrackingMargin, pt); - if (/*!(l->isDead) && */(l->radius <= 2*rmx-5)) {lastActiveLayer = j; last = l; break;} + if (/*!(l->isDead) && */(l->radius <= 2*(rmx-5))) {lastActiveLayer = j; last = l; break;} } if (lastActiveLayer<0) { printf("No active layer with radius < %f is found, pt = %f\n",rmx, pt); return kFALSE; } - // printf("PT=%f 2Rpt=%f Rlr=%f\n",pt,2*rmx,last->radius); // - if (!PropagateToR(&probTr,last->radius + kTrackingMargin,bGauss,1)) return kFALSE; - //if (!probTr.PropagateTo(last->radius,bGauss)) continue; - // reset cov.matrix - // - // rotate to external layer frame - /* - double posL[3]; - probTr.GetXYZ(posL); // lab position - double phiL = TMath::ATan2(posL[1],posL[0]); - if (!probTr.Rotate(phiL)) { - printf("Failed to rotate to the frame (phi:%+.3f)of Extertnal layer at %.2f\n", - phiL,last->radius); - probTr.Print(); - exit(1); + for (int il=1;il<=lastActiveLayer;il++) { + CylLayerK *lr = (CylLayerK*) fLayers.At(il); + AliExternalTrackParam probTrLast(probTr); + bool ok = PropagateToR(&probTrLast,lr->radius,bGauss,1); + if (ok) ok = probTrLast.CorrectForMeanMaterial(lr->radL, 0, mass , kTRUE); + if (ok && lr->xrho>0) { + for (int ise=xrhosteps;ise--;) { + ok = probTrLast.CorrectForMeanMaterial(0, -lr->xrho/xrhosteps, mass , kTRUE); + if (!ok) break; + } } - */ - if (!probTr.Rotate(probTr.Phi())) return kFALSE; // define large errors in track proper frame (snp=0) + if (ok && lr->radius>1e-3 && !lr->isDead) { + ok = probTrLast.Rotate(probTrLast.PhiPos()) && TMath::Abs( probTrLast.GetSnp() )GetName()); Bool_t isVertex = name.Contains("vertex"); + Bool_t isTOF = name.Contains("tof"); // - if (!PropagateToR(&probTr,layer->radius,bGauss,-1)) { - fEfficiency[0] = 0.; - return kFALSE; // exit(1); - } - // if (!probTr.PropagateTo(last->radius,bGauss)) exit(1); // - // rotate to frame with X axis normal to the surface + if (!PropagateToR(&probTr,layer->radius,bGauss,-1)) return kFALSE; //exit(1); if (!isVertex) { double pos[3]; probTr.GetXYZ(pos); // lab position @@ -1556,10 +1564,8 @@ Bool_t DetectorK::SolveTrack(TrackSol& ts) { if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi); if (!probTr.Rotate(phi)) { printf("Failed to rotate to the frame (phi:%+.3f)of layer at %.2f at XYZ: %+.3f %+.3f %+.3f (pt=%+.3f)\n", - phi,layer->radius,pos[0],pos[1],pos[2],pt); - + phi,layer->radius,pos[0],pos[1],pos[2],pt); probTr.Print(); - fEfficiency[0] = 0.; return kFALSE; // exit(1); } } @@ -1569,7 +1575,7 @@ Bool_t DetectorK::SolveTrack(TrackSol& ts) { printf("SaveInw %d (%f) ",j,layer->radius); probTr.Print(); } // - if (!isVertex && !layer->isDead) { + if (!isVertex && !isTOF && !layer->isDead) { // // create fake measurement with the errors assigned to the layer // account for the measurement there @@ -1580,7 +1586,6 @@ Bool_t DetectorK::SolveTrack(TrackSol& ts) { printf("Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}\n", meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]); probTr.Print(); - fEfficiency[0] = 0.; return kFALSE; // exit(1); } } @@ -1589,9 +1594,17 @@ Bool_t DetectorK::SolveTrack(TrackSol& ts) { if (layer->radL>0 && !probTr.CorrectForMeanMaterial(layer->radL, 0, mass , kTRUE)) { printf("Failed to apply material correction, X/X0=%.4f\n",layer->radL); probTr.Print(); - fEfficiency[0] = 0.; return kFALSE; // exit(1); } + if (layer->xrho>0) { // correct in small steps + for (int ise=xrhosteps;ise--;) { + if (!probTr.CorrectForMeanMaterial(0, layer->xrho/xrhosteps, mass , kTRUE)) { + printf("Failed to apply material correction, xrho=%.4f\n",layer->xrho); + probTr.Print(); + return kFALSE; // exit(1); + } + } + } } // // BACKWORD TRACKING +++++++++++++++++ @@ -1606,7 +1619,6 @@ Bool_t DetectorK::SolveTrack(TrackSol& ts) { Bool_t doLikeAliRoot = 0; // don't do the "combined info" but do like in Aliroot - // RESET Covariance Matrix ( to 10 x the estimate -> as it is done in AliExternalTrackParam) // mIstar.UnitMatrix(); // start with unity if (doLikeAliRoot) { @@ -1629,12 +1641,13 @@ Bool_t DetectorK::SolveTrack(TrackSol& ts) { } } //probTr.Rotate(0); - for (Int_t j=0; j<=lastActiveLayer; j++) { // Layer loop + for (Int_t j=0; j<=lastReachedLayer; j++) { // Layer loop // layer = (CylLayerK*)fLayers.At(j); TString name(layer->GetName()); Bool_t isVertex = name.Contains("vertex"); - if (!PropagateToR(&probTr, layer->radius,bGauss,1)) exit(1); + Bool_t isTOF = name.Contains("tof"); + if (!PropagateToR(&probTr, layer->radius,bGauss,1)) return kFALSE;//exit(1); // if (!isVertex) { // rotate to frame with X axis normal to the surface @@ -1663,7 +1676,7 @@ Bool_t DetectorK::SolveTrack(TrackSol& ts) { covCmb[1] = 0; // create fake measurement with the errors assigned to the layer // account for the measurement there - if (!isVertex && !layer->isDead) { + if (!isVertex && !isTOF && !layer->isDead) { double meas[2] = {probTr.GetY(),probTr.GetZ()}; double measErr2[3] = {layer->phiRes*layer->phiRes,0,layer->zRes*layer->zRes}; // @@ -1680,36 +1693,28 @@ Bool_t DetectorK::SolveTrack(TrackSol& ts) { probTr.Print(); return kFALSE; // exit(1); } + if (layer->xrho>0) { // correct in small steps + for (int ise=xrhosteps;ise--;) { + if (!probTr.CorrectForMeanMaterial(0, -layer->xrho/xrhosteps, mass , kTRUE)) { + printf("Failed to apply material correction, xrho=%.4f\n",-layer->xrho); + probTr.Print(); + return kFALSE; // exit(1); + } + } + } // save outward parameters at this layer: after the update new( saveParOutwardA[j] ) AliExternalTrackParam(probTr); // - } - - // here we calculate the probability of adding good hits - for (Int_t j=0; j<=lastActiveLayer; j++) { // Layer loop - // - layer = (CylLayerK*)fLayers.At(j); - TString name(layer->GetName()); - Bool_t isVertex = name.Contains("vertex"); - if (!isVertex && !layer->isDead) { - - // combine covariance matrices before update - Double_t rphiErrorOut = TMath::Sqrt(((AliExternalTrackParam*)saveParOutwardB[j])->GetSigmaY2()); - Double_t rphiErrorIn = TMath::Sqrt(((AliExternalTrackParam*)saveParInward[j])->GetSigmaY2()); - Double_t rphiErrorComb = rphiErrorOut * rphiErrorIn / TMath::Sqrt(rphiErrorOut * rphiErrorOut + rphiErrorIn * rphiErrorIn); - Double_t zErrorOut = TMath::Sqrt(((AliExternalTrackParam*)saveParOutwardB[j])->GetSigmaZ2()); - Double_t zErrorIn = TMath::Sqrt(((AliExternalTrackParam*)saveParInward[j])->GetSigmaZ2()); - Double_t zErrorComb = zErrorOut * zErrorIn / TMath::Sqrt(zErrorOut * zErrorOut + zErrorIn * zErrorIn); - - // calculate probability of adding a good hit - Double_t rphiError = TMath::Sqrt( rphiErrorComb * rphiErrorComb + layer->phiRes * layer->phiRes ); - Double_t zError = TMath::Sqrt( zErrorComb * zErrorComb + layer->zRes * layer->zRes ); - // printf(" --- the efficiency at layer %f --> %f \n", layer->radius, ProbGoodHit( layer->radius, rphiError, zError )); - fEfficiency[0] *= ProbGoodHit( layer->radius, rphiError, zError ); - fLayerEfficiency[j] = ProbGoodHit( layer->radius, rphiError, zError ); + // good hit probability calculation + if (!isVertex && !layer->isDead) { + AliExternalTrackParam* trCmb = (AliExternalTrackParam*)ts.fTrackCmb[j]; + double sigYCmb = TMath::Sqrt(trCmb->GetSigmaY2()+layer->phiRes*layer->phiRes); + double sigZCmb = TMath::Sqrt(trCmb->GetSigmaZ2()+layer->zRes*layer->zRes); + fGoodHitProb[j] = ProbGoodChiSqHit(layer->radius * 100., sigYCmb * 100., sigZCmb * 100.); + if (!isTOF) + fGoodHitProb[0] *= fGoodHitProb[j]; } } - // probTr.SetUseLogTermMS(kFALSE); // Reset of MS term usage to avoid problems since its static // @@ -1764,11 +1769,14 @@ Bool_t DetectorK::CalcITSEff(TrackSol& ts, Bool_t verbose) // if (verbose) { const double kCnv=1e4; - printf("%s:\t%5.1f %.4f %7.0f | %6.0f %6.0f -> %.3f | %6.0f %6.0f -> %.3f | %6.0f %6.0f -> %.3f\n", + printf("%s:\t%5.1f %.4f %7.0f | %6.0f %6.0f -> %.3f | %6.0f %6.0f -> %.3f | %6.0f %6.0f -> %.3f --> %.3f --> %.3f \n", l->GetName(),l->radius,l->radL,HitDensity(l->radius), sigYInw*kCnv,sigZInw*kCnv,probLayInw(2,nITSAct), sigYOut*kCnv,sigZOut*kCnv,probLayOut(2,nITSAct), - sigYCmb*kCnv,sigZCmb*kCnv,probLayCmb(2,nITSAct)); + sigYCmb*kCnv,sigZCmb*kCnv,probLayCmb(2,nITSAct), + ProbGoodHit(l->radius, sigYCmb, sigZCmb), + ProbGoodChiSqHit(l->radius, sigYCmb, sigZCmb) + ); } nITSAct++; ilr++; @@ -2353,6 +2361,15 @@ void DetectorK::MakeStandardPlots(Bool_t add, Int_t color, Int_t linewidth, cons pointResZ->SetName(Form("pointZRes%dadd",0)); pointResZ->Draw("L"); } + if (outGr) { + HistoManager hm("",outGr); + hm.AddGraph(eff); + hm.AddGraph(momRes); + hm.AddGraph(pointResR); + hm.AddGraph(pointResZ); + hm.Write(); + hm.Clear(); + } } @@ -2581,50 +2598,53 @@ Bool_t DetectorK::PropagateToR(AliExternalTrackParam* trc, double r, double b, i double rr = r*r; int iter = 0; const double kTiny = 1e-6; - const Double_t kEpsilon = 0.00001; + const Double_t kEpsilonX = 0.00001, kEpsilonR = 0.01; // if (verboseR) { printf("Prop to %f d=%d ",r,dir); trc->Print(); } - if (!GetXatLabR(trc, r ,xToGo, b, dir)) { - printf("Track with pt=%f cannot reach radius %f\n",trc->Pt(),r); - return kFALSE; - } - - Double_t xpos = trc->GetX(); - dir = (xpos kEpsilon) { - Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep); - Double_t x = xpos+step; - Double_t xyz0[3],xyz1[3],param[7]; - trc->GetXYZ(xyz0); //starting global position - if (!trc->PropagateTo(x,b)) return kFALSE; - xpos = trc->GetX(); - } - // - double rreal = TMath::Sqrt(xpos*xpos+trc->GetY()*trc->GetY()); - // printf("Rtgt=%f Rreal=%f\n",r,rreal); - if (false && r>0.5) { - if (!trc->Rotate(trc->PhiPos())) { - printf("Failed to rotate to layer local frame %f | ",trc->PhiPos()); trc->Print(); + while(1) { + + if (!GetXatLabR(trc, r ,xToGo, b, dir)) { + printf("Track with pt=%f cannot reach radius %f\n",trc->Pt(),r); return kFALSE; } - } - else { - if (!trc->Rotate(trc->Phi())) { - printf("Failed to rotate to track local frame %f | ",trc->Phi()); trc->Print(); - return kFALSE; + + Double_t xpos = trc->GetX(); + dir = (xpos kEpsilonX) { + Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep); + Double_t x = xpos+step; + // Double_t xyz0[3],xyz1[3],param[7]; + // trc->GetXYZ(xyz0); //starting global position + if (!trc->PropagateTo(x,b)) return kFALSE; + xpos = trc->GetX(); + } + // + double drreal = r - TMath::Sqrt(xpos*xpos+trc->GetY()*trc->GetY()); + if (!iter && ((dir>0 && drreal>kEpsilonR) || (dir<0 && drreal<-kEpsilonR)) ) { // apparently the phase changes by more than pi/2 + iter++; + if (!trc->Rotate(trc->Phi())) { + printf("Failed to rotate to track local frame %f in the large phase change mode| ",trc->Phi()); trc->Print(); + return kFALSE; + } + continue; // another iteration + } + // printf("Rtgt=%f Rreal=%f\n",r,rreal); + if (r>0.5) { + if (!trc->Rotate(trc->PhiPos())) { + printf("Failed to rotate to layer local frame %f | ",trc->PhiPos()); trc->Print(); + return kFALSE; + } } + else { + if (!trc->Rotate(trc->Phi())) { + printf("Failed to rotate to track local frame %f | ",trc->Phi()); trc->Print(); + return kFALSE; + } + } + break; } - - // make sure we reached the requested radius - // if not, do it again - double pos[3]; - trc->GetXYZ(pos); // lab position - if (fabs(std::hypot(pos[0], pos[1]) - r) > 0.1) { - return PropagateToR(trc, r, b, dir, maxStep); - } - return kTRUE; } diff --git a/src/DetectorK/DetectorK.h b/src/DetectorK/DetectorK.h index 49d84a2..e262c0a 100644 --- a/src/DetectorK/DetectorK.h +++ b/src/DetectorK/DetectorK.h @@ -6,6 +6,7 @@ #include #include #include +#include "HistoManager.h" /*********************************************************** @@ -113,7 +114,7 @@ class DetectorK : public TNamed { public: DetectorK(); - DetectorK(const char *name, const char *title); + DetectorK(char *name,char *title); virtual ~DetectorK(); enum {kNptBins = 50}; // less then 400 !! @@ -133,9 +134,7 @@ class DetectorK : public TNamed { Float_t GetRadiationLength(char *name); Float_t GetResolution(char *name, Int_t axis=0); Float_t GetLayerEfficiency(char *name); - Float_t GetEfficiency() { return fEfficiency[0]; }; - Float_t GetLayerEfficiency(int i) { return fLayerEfficiency[i]; }; - + void PrintLayout(Bool_t full = kFALSE); void PlotLayout(Int_t plotDead = kTRUE); @@ -231,6 +230,8 @@ class DetectorK : public TNamed { Bool_t IsITSLayer(const TString& lname); + Double_t GetGoodHitProb(Int_t i) { return fGoodHitProb[i]; }; + static Bool_t verboseR; protected: @@ -266,9 +267,9 @@ class DetectorK : public TNamed { Double_t fDetPointRes[kMaxNumberOfDetectors][kNptBins]; // array of rphi resolution per layer Double_t fDetPointZRes[kMaxNumberOfDetectors][kNptBins]; // array of z resolution per layer Double_t fEfficiency[kNptBins]; // efficiency - Double_t fLayerEfficiency[kMaxNumberOfDetectors]; // layer efficiency Double_t fFake[kNptBins]; // fake prob - + Double_t fGoodHitProb[kMaxNumberOfDetectors]; // array of good hit probability per layer + Int_t kDetLayer; // layer for which a few more details are extracted Double_t fResolutionRPhiLay[kNptBins]; // array of rphi resolution Double_t fResolutionZLay[kNptBins]; // array of z resolution diff --git a/src/DetectorK/HistoManager.cxx b/src/DetectorK/HistoManager.cxx new file mode 100644 index 0000000..0fb21a8 --- /dev/null +++ b/src/DetectorK/HistoManager.cxx @@ -0,0 +1,383 @@ +#include "HistoManager.h" +#include "TROOT.h" +#include "TSystem.h" + +ClassImp(HistoManager) + +//_______________________________________________________________ +HistoManager::HistoManager(const char* dirname,const char* fname,Bool_t LOAD,const char* prefix) +{ + // + fNHistos = 0; + fDirName = dirname; + SetFileName(fname); + SetFile(0); + if (LOAD && !fDefName.IsNull()) { + int nh = Load(fname,dirname); + Printf("HistoManager::Load was requested: got %d histos from %s/%s",nh,fname,dirname); + if (prefix && prefix[0]!=0) AddPrefix(prefix); + } + // +} + +//_______________________________________________________________ +HistoManager::~HistoManager() +{ + Delete(); +} + +//_______________________________________________________________ +HistoManager* HistoManager::CreateClone(const char* prefix) +{ + // + HistoManager* hm = (HistoManager*)this->Clone(); + hm->AddPrefix(prefix); + TH1* histo; + for (int i=0;iInheritsFrom("TH1")) { + ((TH1*)obj)->SetDirectory(0); + } + } + hm->fNHistos = fNHistos; + hm->SetFileName(fDefName.Data()); + hm->SetDirName(fDirName.Data()); + return hm; + // +} + +//_______________________________________________________________ +Int_t HistoManager::AddHisto(TH1* histo,Int_t at) +{ + // Add new histo either to next free slot (at<0) or to requested position + if (at<0) at = fNHistos; + AddAtAndExpand(histo,at); + histo->SetDirectory(0); + histo->SetUniqueID(at+1); + return fNHistos++; + // +} + +//_______________________________________________________________ +Int_t HistoManager::AddGraph(TGraph* gr,Int_t at) +{ + // Add new histo either to next free slot (at<0) or to requested position + if (at<0) at = fNHistos; + AddAtAndExpand(gr,at); + //histo->SetDirectory(0); + gr->SetUniqueID(at+1); + return fNHistos++; + // +} + +//_______________________________________________________________ +void HistoManager::Compress() +{ + TObjArray::Compress(); + TObject* histo; + for (int i=0;iSetUniqueID(i+1); +} + +//_______________________________________________________________ +void HistoManager::Write(TFile* file) +{ + // Write all histograms to file + if (!fNHistos) return; + Bool_t localfile = kFALSE; + TH1* histo=0; + TFile *lfile=0; + const char* str=0; + if (file) lfile = file; + else { + // Check if the file is not already open + TFile *tmpF = (TFile*)gROOT->GetListOfFiles()->FindObject(fDefName.Data()); + if (tmpF && tmpF->IsOpen()) { + TString opt = tmpF->GetOption(); opt.ToLower(); + if (!opt.Contains("read")) { + lfile = tmpF; + tmpF->cd(); + } + } + } + TString pwd = gDirectory->GetPath(); + if (!lfile) { // have to open + str = fDefName.Data(); + if (!str || !str[0] || str[0] == ' ') fDefName = "histoman"; + if (!fDefName.Contains(".root")) fDefName += ".root"; + lfile = TFile::Open(fDefName.Data(),"UPDATE"); + fDefName = str; + localfile = kTRUE; + } + // + lfile->cd(); + // Create directory (if necessary) + str = fDirName.Data(); + if (str && str[0] && str[0] != ' ') { + if (!lfile->Get(str)) lfile->mkdir(str); + lfile->cd(str); + } + Printf("Writing histogrames to: %s%s",lfile->GetPath(),str); + for (int i=0;i(obj); + TDirectory* dr = 0; + if (histo) { + dr = histo->GetDirectory(); + histo->SetDirectory(0); + } + obj->Write(0,TObject::kOverwrite); + if (dr && histo) histo->SetDirectory(dr); + } + if (localfile) {lfile->Close(); delete lfile;} + TDirectory* oldDir = ((TDirectory *)gROOT->GetListOfFiles()->FindObject(pwd.Data())); + if (oldDir) oldDir->cd(); + // +} + +//_______________________________________________________________ +void HistoManager::Clear(Option_t*) +{ + int nent = GetLast()+1; + for (int i=0;iPrint(option); + // + } + Printf("\nTotal number of defined Histograms: %d",fNHistos); + Printf("\nCurrent output path: %s/%s",fDefName.Data(),fDirName.Data()); +} + +//_______________________________________________________________ +void HistoManager::AddPrefix(const char* pref) +{ + TString prfs = pref; + if (prfs.IsNull()) return; + int nent = GetLast()+1; + for (int i=0;iGetName(); + if (hh->InheritsFrom("TNamed")) ((TNamed*)hh)->SetName(prfs.Data()); + // prfs = pref; + // prfs += hh->GetTitle(); + // hh->SetTitle(prfs.Data()); + } +} + +//_______________________________________________________________ +void HistoManager::AddHistos(HistoManager* hm,Double_t c1) +{ + int nent = GetLast()+1; + int nent1 = hm->GetLast()+1; + if (nent!=nent1) {Error("AddHistos","HistoManagers have different content: %d vs %d",nent,nent1);return;} + for (int i=0;iGetHisto(i); + if (!hh1 || !hh2) continue; + hh1->Add(hh2,c1); + } +} + +//_______________________________________________________________ +void HistoManager::DivideHistos(HistoManager* hm) +{ + int nent = GetLast()+1; + int nent1 = hm->GetLast()+1; + if (nent!=nent1) {Error("DivideHistos","HistoManagers have different content: %d vs %d",nent,nent1);return;} + for (int i=0;iGetHisto(i); + if (!hh1 || !hh2) continue; + hh1->Divide(hh2); + } +} + +//_______________________________________________________________ +void HistoManager::MultiplyHistos(HistoManager* hm) +{ + int nent = GetLast()+1; + int nent1 = hm->GetLast()+1; + if (nent!=nent1) {Error("MultiplyHistos","HistoManagers have different content: %d vs %d",nent,nent1);return;} + for (int i=0;iGetHisto(i); + if (!hh1 || !hh2) continue; + hh1->Multiply(hh2); + } +} + +//_______________________________________________________________ +void HistoManager::ScaleHistos(Double_t c1) +{ + int nent = GetLast()+1; + for (int i=0;iScale(c1); + } +} + +//_______________________________________________________________ +void HistoManager::Sumw2() +{ + int nent = GetLast()+1; + for (int i=0;i(UncheckedAt(i)); + if (hh1) hh1->Sumw2(); + } +} + +//_______________________________________________________________ +void HistoManager::SetFile(TFile* file) +{ + if (file) fDefName = file->GetName(); +} + +//_______________________________________________________________ +void HistoManager::DelHisto(Int_t at) +{ + TH1* hist = GetHisto(at); + if (hist) { + RemoveAt(at); + delete hist; + } +} + +//_______________________________________________________________ +void HistoManager::Purify(Bool_t emptyToo) +{ + // remove empty slots, optionally removing empty histos too + int last = GetLast()+1; + if (emptyToo) + for (int i=0;iGetEntries()<1) { DelHisto(i); fNHistos--;} + } + Compress(); + // +} + +//_____________________________________________________________________________ +void HistoManager::SetFileName(const char *name) +{fDefName = name; gSystem->ExpandPathName(fDefName);} + +//_____________________________________________________________________________ +void HistoManager::Reset() +{ + int last = GetLast()+1; + for (int i=0;iReset(); + } +} + +//_____________________________________________________________________________ +Int_t HistoManager::Load(const char* fname,const char* dirname) +{ + TString flpath = fname; + gSystem->ExpandPathName(flpath); + TFile* file = TFile::Open(flpath.Data()); + if (!file) {Printf("Error: no file %s",fname); return 0;} + if (dirname && dirname[0] && dirname[0] != ' ') { + if (!file->Get(dirname)) { + Printf("Error: no %s directory in file %s",dirname,fname); + file->Close(); delete file; + return 0; + } + else file->cd(dirname); + } + // + int count = 0; + TList* lst = gDirectory->GetListOfKeys(); + TIter next(lst); + TObject* obj; + while ((obj=next())) { + if (FindObject(obj->GetName())) continue; // already added + TObject* hst = gDirectory->Get(obj->GetName()); + int ID = hst->GetUniqueID(); + TH1* h = dynamic_cast(hst); + if (h) { + AddHisto(h, ID-1); + count++; + continue; + } + TGraph* gr = dynamic_cast(hst); + if (gr) { + AddGraph(gr, ID-1); + count++; + continue; + } + } + file->Close(); + delete file; + return count; +} + +//_____________________________________________________________________________ +void HistoManager::SetColor(Color_t tcolor) +{ + int last = GetLast()+1; + for (int i=0;iSetLineColor(tcolor); + hist->SetMarkerColor(tcolor); + } +} + +//_____________________________________________________________________________ +void HistoManager::SetMarkerStyle(Style_t mstyle,Size_t msize) +{ + int last = GetLast()+1; + for (int i=0;iSetMarkerStyle(mstyle); + hist->SetMarkerSize(msize); + } +} + +//_____________________________________________________________________________ +void HistoManager::SetMarkerSize(Size_t msize) +{ + int last = GetLast()+1; + for (int i=0;iSetMarkerSize(msize); + } +} + diff --git a/src/DetectorK/HistoManager.h b/src/DetectorK/HistoManager.h new file mode 100644 index 0000000..45be301 --- /dev/null +++ b/src/DetectorK/HistoManager.h @@ -0,0 +1,74 @@ +#ifndef HISTOMANAGER_H +#define HISTOMANAGER_H + +#include "TH1.h" +#include "TH2.h" +#include "TGraph.h" +#include "TProfile.h" +#include "TFile.h" +#include "TObjArray.h" +class TROOT; +class TSystem; + +class HistoManager +: public TObjArray +{ + //#ifdef USE_USING + using TCollection::Print; + using TCollection::Write; + //#endif + public: + HistoManager(const char* dirname="",const char* fname="histoman.root",Bool_t LOAD=kFALSE,const char* prefix=""); + ~HistoManager(); + HistoManager* CreateClone(const char* prefix); + // + Int_t GetNHistos() const {return fNHistos;} + TGraph* GetGraph(Int_t id) const {return id<=GetLast() ? dynamic_cast(UncheckedAt(id)) : 0;} + TH1* GetHisto(Int_t id) const {return id<=GetLast() ? dynamic_cast(UncheckedAt(id)) : 0;} + TH1* GetHisto(char* name) const {return dynamic_cast( FindObject(name) );} + TH1F* GetHisto1F(Int_t id) const {return dynamic_cast( UncheckedAt(id) );} + TH2F* GetHisto2F(Int_t id) const {return dynamic_cast( UncheckedAt(id) );} + TProfile* GetHistoP(Int_t id) const {return dynamic_cast( UncheckedAt(id) );} + Int_t AddHisto(TH1* histo,Int_t at=-1); + Int_t AddGraph(TGraph* gr,Int_t at=-1); + void DelHisto(Int_t at); + void SetFile(TFile* file); + void SetFileName(const char* fname); + char* GetFileName() const {return (char*) fDefName.Data();} + void SetDirName(const char* name) {fDirName = name;} + char* GetDirName() const {return (char*) fDirName.Data();} + // + void Reset(); + void Write(TFile* file=0); + Int_t Write(const char* flname, int =0, int =0) {SetFileName(flname);Write();return 0;} + void AddPrefix(const char* pref); + void AddHistos(HistoManager* hm, Double_t c1 = 1); + void DivideHistos(HistoManager* hm); + void MultiplyHistos(HistoManager* hm); + void ScaleHistos(Double_t c1 = 1); + void SetColor(Color_t tcolor = 1); + void SetMarkerStyle(Style_t mstyle = 1, Size_t msize = 1); + void SetMarkerSize(Size_t msize = 1); + void Sumw2(); + Int_t Load(const char* fname,const char* dirname=""); + // + void Purify(Bool_t emptyToo=kFALSE); + void Print(Option_t* option="") const; + void Clear(Option_t* option=""); + void Delete(Option_t* option=""); + virtual void Compress(); + // + private: + Int_t fNHistos; // Number of Histogrames defined + TString fDefName; // Default file name + TString fDirName; // Directory name in the output file + // + ClassDef(HistoManager,0) // NA60 Histograms manager +}; + + + + + + +#endif From fe461f81864d4d3d3eafe4289118c7353a78502e Mon Sep 17 00:00:00 2001 From: Roberto Preghenella Date: Sat, 31 Jul 2021 17:06:24 +0200 Subject: [PATCH 09/10] Adjust fwdparam routine to match current code --- src/lutWrite.cc | 34 ++++++++++++++++++++-------------- 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/src/lutWrite.cc b/src/lutWrite.cc index e81746c..90fc7fa 100644 --- a/src/lutWrite.cc +++ b/src/lutWrite.cc @@ -8,6 +8,8 @@ DetectorK fat; void diagonalise(lutEntry_t &lutEntry); static float etaMaxBarrel = 1.75; +bool usePara = true; // use fwd parameterisation + bool fatSolve(lutEntry_t &lutEntry, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, int itof = 0, int otof = 0, int q = 1) { @@ -37,15 +39,19 @@ fwdSolve(float *covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000) } bool -fwdPara(float *covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, float Bfield = 0.5, int layer = 0, int what = 0) +fwdPara(lutEntry_t &lutEntry, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, float Bfield = 0.5) { + lutEntry.valid = false; + // parametrised forward response; interpolates between FAT at eta = 1.75 and a fixed parametrisation at eta = 4; only diagonal elements - float eff = 0; - float covmbarrel[15] = {0}; if (fabs(eta) < etaMaxBarrel || fabs(eta) > 4) return false; - if (!fatSolve(&eff, covmbarrel, pt, etaMaxBarrel, mass, layer, what)) return false; + if (!fatSolve(lutEntry, pt, etaMaxBarrel, mass)) return false; + float covmbarrel[15] = {0}; + for (int i = 0; i < 15; ++i) { + covmbarrel[i] = lutEntry.covm[i]; + } // parametrisation at eta = 4 double beta = 1./sqrt(1+mass*mass/pt/pt/cosh(eta)/cosh(eta)); @@ -76,15 +82,15 @@ fwdPara(float *covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, f // Fill cov matrix diag for (int i = 0; i < 15; ++i) - covm[i] = 0; - - covm[0] = covmbarrel[0]; - if (dcaxy2 > covm[0]) covm[0] = dcaxy2; - covm[2] = covmbarrel[2]; - if (dcaz2 > covm[2]) covm[2] = dcaz2; - covm[5] = covmbarrel[5]; // sigma^2 sin(phi) - covm[9] = covmbarrel[9]; // sigma^2 tanl - covm[14] = momres_tot*momres_tot/pt/pt/pt/pt; // sigma^2 1/pt + lutEntry.covm[i] = 0; + + lutEntry.covm[0] = covmbarrel[0]; + if (dcaxy2 > lutEntry.covm[0]) lutEntry.covm[0] = dcaxy2; + lutEntry.covm[2] = covmbarrel[2]; + if (dcaz2 > lutEntry.covm[2]) lutEntry.covm[2] = dcaz2; + lutEntry.covm[5] = covmbarrel[5]; // sigma^2 sin(phi) + lutEntry.covm[9] = covmbarrel[9]; // sigma^2 tanl + lutEntry.covm[14] = momres_tot*momres_tot/pt/pt/pt/pt; // sigma^2 1/pt return true; } @@ -167,7 +173,7 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, lutEntry.eff = 1.; bool retval = true; if (usePara) { - retval = fwdPara(lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass, field, layer, what); + retval = fwdPara(lutEntry, lutEntry.pt, lutEntry.eta, lutHeader.mass, field); } else { retval = fwdSolve(lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass); From bb9b7986fcc6678f9f9ae9f7e0db948abcc8612b Mon Sep 17 00:00:00 2001 From: Roberto Preghenella Date: Sun, 1 Aug 2021 15:48:35 +0200 Subject: [PATCH 10/10] Add new definition for the efficiency --- examples/smearing/lutRead.C | 2 ++ src/TrackSmearer.cc | 9 +++++++-- src/TrackSmearer.hh | 2 ++ src/lutCovm.hh | 3 ++- src/lutWrite.cc | 23 +++++++++++++++++++++-- 5 files changed, 34 insertions(+), 5 deletions(-) diff --git a/examples/smearing/lutRead.C b/examples/smearing/lutRead.C index d43fee4..9226aa1 100644 --- a/examples/smearing/lutRead.C +++ b/examples/smearing/lutRead.C @@ -1,5 +1,6 @@ enum EWhat { kEfficiency, + kEfficiency2, kEfficiencyInnerTOF, kEfficiencyOuterTOF, kPtResolution, @@ -44,6 +45,7 @@ lutRead(int pdg, const char *filename, int what, int vs, float nch = 0., float r if (vs == kPt) cen = lutEntry->pt; double val = 0.; if (what == kEfficiency) val = lutEntry->eff * 100.; // efficiency (%) + if (what == kEfficiency2) val = lutEntry->eff2 * 100.; // efficiency (%) if (what == kEfficiencyInnerTOF) val = lutEntry->itof * 100.; // efficiency (%) if (what == kEfficiencyOuterTOF) val = lutEntry->otof * 100.; // efficiency (%) if (what == kPtResolution) val = sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.; // pt resolution (%) diff --git a/src/TrackSmearer.cc b/src/TrackSmearer.cc index 9a17fa8..598164d 100644 --- a/src/TrackSmearer.cc +++ b/src/TrackSmearer.cc @@ -96,8 +96,13 @@ bool TrackSmearer::smearTrack(O2Track &o2track, lutEntry_t *lutEntry) { // generate efficiency - if (mUseEfficiency && (gRandom->Uniform() > lutEntry->eff)) - return false; + if (mUseEfficiency) { + auto eff = 0.; + if (mWhatEfficiency == 1) eff = lutEntry->eff; + if (mWhatEfficiency == 2) eff = lutEntry->eff2; + if (gRandom->Uniform() > eff) + return false; + } // transform params vector and smear double params_[5]; for (int i = 0; i < 5; ++i) { diff --git a/src/TrackSmearer.hh b/src/TrackSmearer.hh index b6af3f8..4f3bc70 100644 --- a/src/TrackSmearer.hh +++ b/src/TrackSmearer.hh @@ -25,6 +25,7 @@ public: /** LUT methods **/ bool loadTable(int pdg, const char *filename, bool forceReload = false); void useEfficiency(bool val) { mUseEfficiency = val; }; + void setWhatEfficiency(int val) { mWhatEfficiency = val; }; lutHeader_t *getLUTHeader(int pdg) { return mLUTHeader[getIndexPDG(pdg)]; }; lutEntry_t *getLUTEntry(int pdg, float nch, float radius, float eta, float pt); @@ -53,6 +54,7 @@ protected: lutHeader_t *mLUTHeader[nLUTs] = {nullptr}; lutEntry_t *****mLUTEntry[nLUTs] = {nullptr}; bool mUseEfficiency = true; + int mWhatEfficiency = 1; float mdNdEta = 1600.; }; diff --git a/src/lutCovm.hh b/src/lutCovm.hh index 8cb8c6a..acceefd 100644 --- a/src/lutCovm.hh +++ b/src/lutCovm.hh @@ -2,7 +2,7 @@ /// @email: preghenella@bo.infn.it #pragma once -#define LUTCOVM_VERSION 20210731 +#define LUTCOVM_VERSION 20210801 struct map_t { int nbins = 1; @@ -56,6 +56,7 @@ struct lutEntry_t { float pt = 0.; bool valid = false; float eff = 0.; + float eff2 = 0.; float itof = 0.; float otof = 0.; float covm[15] = {0.}; diff --git a/src/lutWrite.cc b/src/lutWrite.cc index 90fc7fa..371bb1a 100644 --- a/src/lutWrite.cc +++ b/src/lutWrite.cc @@ -23,11 +23,28 @@ fatSolve(lutEntry_t &lutEntry, float pt = 0.1, float eta = 0.0, float mass = 0.1 if (!trPtr) return false; lutEntry.valid = true; - lutEntry.eff = fat.GetGoodHitProb(0); lutEntry.itof = fat.GetGoodHitProb(itof); lutEntry.otof = fat.GetGoodHitProb(otof); for (int i = 0; i < 15; ++i) lutEntry.covm[i] = trPtr->GetCovariance()[i]; - + + // define the efficiency + auto totfake = 0.; + lutEntry.eff = 1.; + for (int i = 1; i < 20; ++i) { + auto igoodhit = fat.GetGoodHitProb(i); + if (igoodhit <= 0. || i == itof || i == otof) continue; + lutEntry.eff *= igoodhit; + auto pairfake = 0.; + for (int j = i + 1; j < 20; ++j) { + auto jgoodhit = fat.GetGoodHitProb(j); + if (jgoodhit <= 0. || j == itof || j == otof) continue; + pairfake = (1. - igoodhit) * (1. - jgoodhit); + break; + } + totfake += pairfake; + } + lutEntry.eff2 = (1. - totfake); + return true; } @@ -164,6 +181,7 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, // printf(" --- fatSolve: error \n"); lutEntry.valid = false; lutEntry.eff = 0.; + lutEntry.eff2 = 0.; for (int i = 0; i < 15; ++i) lutEntry.covm[i] = 0.; } @@ -171,6 +189,7 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2, else { // printf(" --- fwdSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); lutEntry.eff = 1.; + lutEntry.eff2 = 1.; bool retval = true; if (usePara) { retval = fwdPara(lutEntry, lutEntry.pt, lutEntry.eta, lutHeader.mass, field);