diff --git a/StDb/idl/pidCorrection.idl b/StDb/idl/pidCorrection.idl new file mode 100644 index 00000000000..5dfd9de8831 --- /dev/null +++ b/StDb/idl/pidCorrection.idl @@ -0,0 +1,26 @@ +//:Description: Drift Distance depended correction +//:Synonyms:::: +//:Source: +//:Update: +//:Update frequncy: +//:Reminder: +//:Recall frequency: +//:Size of Data: +//:Pointer to data: pidCorrection.time: +struct pidCorrection { + long idx; // row index > 0 if it is real + long nrows; // total no. of real rows in the table; For Db interface (where nrows = 50) + long type; // type = 0 polymonical fit, use only [min,max] + long var; // fit variable: 0 => pmomL10, 1 => bgL10, + long particle; // from StEvent/StPidParticleDefinition.h : kUndef = -1, kPidElectron = 0, Proton = 1, Kaon = 2, Pion = 3, Muon = 4, Deuteron = 5, Triton = 6, + // He3 = 7, Alpha = 8, He6 = 9, Li5 = 10, Li6,= 11, Li7 = 12, Be7 = 13, Be9 = 14, Be10 = 15, B11 = 16 + long charge; // +/-1, 0 does not matter + long pull; // != 0 calculated pull correction, == 0 to value + long det; // from StdEdxY2Maker/StTrackPiD.h : kUndef = 0, kI70 = 1, kI70U = 2, kFit = 3, kFitU = 4, kdNdx = 5, kdNdxU = 6, kBTof -7 , kETof = 8, kMtd = 9, kBEmc = 10 + long npar; // npar < 0, X = exp(x) paramterization; abs(npar) >= 100 cut on range [min.max] + double OffSet; // for future use + double min; // fit range + double max; // + double a[10]; // a[npar] + char comment[32]; +}; diff --git a/StRoot/StBichsel/StBichselLinkDef.h b/StRoot/StBichsel/StBichselLinkDef.h index ed15be8b35f..de843e256e3 100644 --- a/StRoot/StBichsel/StBichselLinkDef.h +++ b/StRoot/StBichsel/StBichselLinkDef.h @@ -12,4 +12,7 @@ #pragma link C++ function StdEdxPull::EvalPred2(Double_t, Double_t, UChar_t, Int_t); #pragma link C++ function StdEdxPull::EvalDeV2(Double_t, Double_t, Double_t, UChar_t, Int_t); #pragma link C++ function StdEdxPull::Eval2(Double_t, Double_t, Double_t, UChar_t, Int_t); + +#pragma link C++ class St_spline3-; +#pragma link C++ class spline3_st+; #endif diff --git a/StRoot/StBichsel/St_spline3C.cxx b/StRoot/StBichsel/St_spline3C.cxx new file mode 100644 index 00000000000..9c891391a6c --- /dev/null +++ b/StRoot/StBichsel/St_spline3C.cxx @@ -0,0 +1,73 @@ +#include "Riostream.h" +#include "St_spline3C.h" +#include "TString.h" +#include "TInterpreter.h" +#include "TSystem.h" +//________________________________________________________________________________ +St_spline3 *St_spline3C::Open(const Char_t *path) { + St_spline3 *table = 0; + TString PATH(path); + TString Dir(gSystem->DirName(PATH)); + TString File(gSystem->BaseName(PATH)); + File += ".C"; + TString pathF(".:./StarDb/"); pathF += Dir + ":$STAR/StarDb/" + Dir; + Char_t *file = gSystem->Which(pathF,File,kReadPermission); + if (! file) { + std::cout << Form("Fatal::St_spline3C::Open \tFile %s has not been found in path %s",File.Data(),pathF.Data()) << std::endl; + return table; + } else { + std::cout << Form("Warning::St_spline3C::Open \tFile %s has been found as %s",File.Data(),file) << std::endl; + } + TString command(".L "); command += file; TInterpreter::EErrorCode ee; + gInterpreter->ProcessLine(command,&ee); + if (ee) { //assert(!ee); + std::cout << Form("Fatal::St_spline3C::Open has failed to read \tFile %s",file) << std::endl; + delete [] file; + return table; + } + table = (St_spline3 *) gInterpreter->Calc("CreateTable()",&ee); + if (! table) {//assert(table); + std::cout << Form("Fatal::St_spline3C::Open has failed to load \tFile %s",file) << std::endl; + delete [] file; + return table; + } + table->Print(0,1); + command.ReplaceAll(".L ",".U "); + gInterpreter->ProcessLine(command,&ee); + if (ee) { // assert(!ee); + std::cout << Form("Fatal::St_spline3C::Open has failed to unload \tFile %s",file) << std::endl; + delete [] file; + SafeDelete(table); + return table; + } + return table; +} +//________________________________________________________________________________ +St_spline3C::St_spline3C(St_spline3 *table) : TChair(table), fSpline(0), fFunc(0), fValid(kTRUE) { + if (table) { + fSpline = new TSpline3("Spline3", Xknots(), Yknots(), nknots(), option(), ValBeg(), ValEnd()); + fSpline->SetLineColor(2); + fXmin = Xknots()[0] - 0.1; + fXmax = Xknots()[nknots()-1] + 0.1; + fFunc = new TF1(GetName(), this, fXmin, fXmax, 0, "St_spline3C"); + fFunc->SetNpx(100); + fFunc->Save(Xknots()[0], Xknots()[nknots()-1], 0., 0., 0., 0.); + } else { + fValid = kFALSE; + } +} +#define MakeChairInstance3(CLASS,PATH) \ + ClassImp(CLASS); \ + CLASS *CLASS::fgInstance = 0; \ + CLASS *CLASS::instance() { \ + if (fgInstance && ! fgInstance->IsValid()) return 0; \ + if (fgInstance) return fgInstance; \ + St_spline3 *table = St_spline3C::Open(# PATH); \ + fgInstance = new CLASS(table); \ + return fgInstance; \ + } +//________________________________________________________________________________ +MakeChairInstance3(Stspline3LndNdxL10,dEdxModel/spline3LndNdxL10); +MakeChairInstance3(StElectonsDEV_dEdx,dEdxModel/ElectonsDEV_dEdx); +MakeChairInstance3(StPionDEV_dEdx,dEdxModel/PionDEV_dEdx); +MakeChairInstance3(StProtonDEV_dEdx,dEdxModel/ProtonDEV_dEdx); diff --git a/StRoot/StBichsel/St_spline3C.h b/StRoot/StBichsel/St_spline3C.h new file mode 100644 index 00000000000..c62ee923f07 --- /dev/null +++ b/StRoot/StBichsel/St_spline3C.h @@ -0,0 +1,77 @@ +#ifndef St_spline3C_h +#define St_spline3C_h + +#include "TChair.h" +//#include "tables/St_spline3_Table.h" +#include "St_spline3_Table.h" +#include "TSpline.h" +#include "TF1.h" +//________________________________________________________________________________ +class St_spline3C : public TChair { + public: + St_spline3C(St_spline3 *table=0); + virtual ~St_spline3C() {} + spline3_st *Struct(Int_t i = 0) const {return ((St_spline3*) Table())->GetTable()+i;} + UInt_t getNumRows() const {return GetNRows();} + Int_t nknots(Int_t i = 0) const {return Struct(i)->nknots;} + Double_t* Xknots(Int_t i = 0) const {return Struct(i)->Xknots;} + Double_t* Yknots(Int_t i = 0) const {return Struct(i)->Yknots;} + Double_t ValBeg(Int_t i = 0) const {return Struct(i)->ValBeg;} + Double_t ValEnd(Int_t i = 0) const {return Struct(i)->ValEnd;} + Char_t* option(Int_t i = 0) const {return (Char_t *) Struct(i)->option;} + TSpline3* Spline() {return fSpline;} + Double_t operator() (Double_t *x, Double_t *p) const {return fSpline->Eval(x[0]);} + TF1* Func() {return fFunc;} + Bool_t IsValid() {return fValid;} + static St_spline3 *Open(const Char_t *path); + Bool_t InRange(Double_t x) {return fXmin <= x && x <= fXmax;} + private: + TSpline3* fSpline; + TF1* fFunc; + Bool_t fValid; + Double_t fXmin; + Double_t fXmax; + ClassDefChair(St_spline3, spline3_st ) + ClassDef(St_spline3C,1) //C++ TChair for spline3 table class +}; +//________________________________________________________________________________ +class Stspline3LndNdxL10 : public St_spline3C {// Log(dN/dx) versus log10(beta*gamma) + public: + static Stspline3LndNdxL10* instance(); + Stspline3LndNdxL10(St_spline3 *table=0) : St_spline3C(table) {} + virtual ~Stspline3LndNdxL10() {fgInstance = 0;} + private: + static Stspline3LndNdxL10* fgInstance; + ClassDef(Stspline3LndNdxL10,1) //C++ TChair for spline3LndNdxL10 +}; +//________________________________________________________________________________ +class StElectonsDEV_dEdx : public St_spline3C {// Log(dN/dx) versus log10(beta*gamma) + public: + static StElectonsDEV_dEdx* instance(); + StElectonsDEV_dEdx(St_spline3 *table=0) : St_spline3C(table) {} + virtual ~StElectonsDEV_dEdx() {fgInstance = 0;} + private: + static StElectonsDEV_dEdx* fgInstance; + ClassDef(StElectonsDEV_dEdx,1) //C++ TChair for ElectonsDEV_dEdx +}; +//________________________________________________________________________________ +class StPionDEV_dEdx : public St_spline3C {// Log(dN/dx) versus log10(beta*gamma) + public: + static StPionDEV_dEdx* instance(); + StPionDEV_dEdx(St_spline3 *table=0) : St_spline3C(table) {} + virtual ~StPionDEV_dEdx() {fgInstance = 0;} + private: + static StPionDEV_dEdx* fgInstance; + ClassDef(StPionDEV_dEdx,1) //C++ TChair for PionDEV_dEdx +}; +//________________________________________________________________________________ +class StProtonDEV_dEdx : public St_spline3C {// Log(dN/dx) versus log10(beta*gamma) + public: + static StProtonDEV_dEdx* instance(); + StProtonDEV_dEdx(St_spline3 *table=0) : St_spline3C(table) {} + virtual ~StProtonDEV_dEdx() {fgInstance = 0;} + private: + static StProtonDEV_dEdx* fgInstance; + ClassDef(StProtonDEV_dEdx,1) //C++ TChair for ProtonDEV_dEdx +}; +#endif diff --git a/StRoot/StBichsel/St_spline3_Table.cxx b/StRoot/StBichsel/St_spline3_Table.cxx new file mode 100644 index 00000000000..c3223c6da91 --- /dev/null +++ b/StRoot/StBichsel/St_spline3_Table.cxx @@ -0,0 +1,11 @@ + +#include "St_spline3_Table.h" +///////////////////////////////////////////////////////////////////////// +// +// Class St_spline3 wraps the STAF table spline3 +// It has been generated "by automatic". Please don't change it "by hand" +// +///////////////////////////////////////////////////////////////////////// + +#include "Stypes.h" +TableImpl(spline3) diff --git a/StRoot/StBichsel/St_spline3_Table.h b/StRoot/StBichsel/St_spline3_Table.h new file mode 100644 index 00000000000..6f60ccef010 --- /dev/null +++ b/StRoot/StBichsel/St_spline3_Table.h @@ -0,0 +1,23 @@ + +#ifndef STAF_St_spline3_Table +#define STAF_St_spline3_Table + +#include "TTable.h" + +#include "spline3.h" + +/*! + * \class St_spline3 + * \brief C++ wrapper for StAF table + * \author Automatic Generation + * \date Fri May 26 18:56:01 2023 + * + * This was generated for version '.DEV2' + */ +class St_spline3 : public TTable +{ + public: + ClassDefTable(St_spline3,spline3_st) + ClassDef(St_spline3,2) //C++ wrapper for StAF table +}; +#endif diff --git a/StRoot/StBichsel/StdEdxModel.cxx b/StRoot/StBichsel/StdEdxModel.cxx index 781498ee046..9022923c133 100644 --- a/StRoot/StBichsel/StdEdxModel.cxx +++ b/StRoot/StBichsel/StdEdxModel.cxx @@ -11,6 +11,7 @@ #include "StdEdxModel.h" #include "Bichsel.h" #include "TCallf77.h" +#include "St_spline3C.h" #ifndef WIN32 #define ggderiv ggderiv_ #else @@ -30,46 +31,62 @@ StdEdxModel* StdEdxModel::instance() { return fgStdEdxModel; } //________________________________________________________________________________ -StdEdxModel::StdEdxModel() : mdNdx(0), fScale(1) +StdEdxModel::StdEdxModel() : fScale(1) , fTmaxL10eV(5) // Tcut = 100 keV - , fGGaus(0), fGausExp(0) - , fpol2F(0), fpol5F(0), fpol6F(0), fpol7F(0) - , fLogkeVperElectron(0) { // LOG_INFO << "StdEdxModel:: use StTpcRSMaker model for dE/dx calculations" << endm; cout << "StdEdxModel:: use StTpcRSMaker model for dE/dx calculations" << endl; - if (! fgStdEdxModel) { - TDirectory *dir = gDirectory; - fgStdEdxModel = this; - const Char_t *path = ".:./StarDb/dEdxModel:$STAR/StarDb/dEdxModel"; - const Char_t *Files[1] = {"dNdx_Bichsel.root"}; - for (Int_t i = 0; i < 1; i++) { // files - Char_t *file = gSystem->Which(path,Files[i],kReadPermission); - if (! file) Fatal("StdEdxModel","File %s has not been found in path %s",Files[i],path); - else Warning("StdEdxModel","File %s has been found as %s",Files[i],file); - TFile *pFile = new TFile(file); - mdNdx = (TH1D *) pFile->Get("dNdx"); if (mdNdx) mdNdx->SetDirectory(0); - assert(mdNdx); - delete pFile; - delete [] file; + memset(beg, 0, end-beg+1); + TDirectory *dir = gDirectory; + fgStdEdxModel = this; + Int_t NFiles = 2; // old + if ( ! Stspline3LndNdxL10::instance()) NFiles = 1; + if ( Stspline3LndNdxL10::instance()->IsValid()) NFiles = 1; + + const Char_t *path = ".:./StarDb/dEdxModel:$STAR/StarDb/dEdxModel"; +#ifndef __HEED__ + const Char_t *Files[2] = {"dNdE_Bichsel.root","dNdx_Bichsel.root"}; +#else + const Char_t *Files[2] = {"dNdx_Heed.root","dNdx_Heed.root"}; +#endif +#define __Warn_Hist__(__HIST__) {m ## __HIST__ = (TH1D *) pFile->Get(#__HIST__); \ + if (m ## __HIST__) {Warning("StdEdxModel","Histogram %s/%s has been found",m ## __HIST__->GetName(),m ## __HIST__->GetTitle()); m ## __HIST__->SetDirectory(0);}} + for (Int_t i = 0; i < NFiles; i++) { // files + Char_t *file = gSystem->Which(path,Files[i],kReadPermission); + if (! file) Fatal("StdEdxModel","File %s has not been found in path %s",Files[i],path); + else Warning("StdEdxModel","File %s has been found as %s",Files[i],file); + TFile *pFile = new TFile(file); + if (i == 1) { + __Warn_Hist__(dNdx); + __Warn_Hist__(dNdxL10); + __Warn_Hist__(LndNdxL10); + __Warn_Hist__(LndNdxL10Smooth); + mLndNdxL10Spline5 = (TSpline5 *) pFile->Get("LndNdxL10Spline5"); if (mLndNdxL10Spline5) {Warning("StdEdxModel","TSpline5 %s has been found",mLndNdxL10Smooth->GetName());} + assert(mdNdx || mdNdxL10); + } else { + mdNdEL10 = (TH1D *) pFile->Get("dNdEL10"); assert(mdNdEL10); Warning("StdEdxModel","Histogram %s/%s has been found",mdNdEL10->GetName(),mdNdEL10->GetTitle()); mdNdEL10->SetDirectory(0); } - dir->cd(); + delete pFile; + delete [] file; } +#undef __Warn_Hist__ + dir->cd(); fGGaus = new TF1("GGaus",ggaus, -1., 5., 5); fGGaus->SetParNames("NormL","mu","sigma","alpha","k"); fGGaus->SetParameters(0,0,0.3,0.1,0); fGGaus->FixParameter(4,0.0); - + fGausExp = new TF1("GausExp",gausexp, -5., 5., 5); fGausExp->SetParNames("NormL","mu","sigma","k","l"); fGausExp->SetParameters(0,0,0.3,5.,0); fGausExp->FixParameter(4,0.0); InitPar(); // Set normalization point the same as for I70 (increase energy per conduction electron from 20 eB to 52 eV) - Double_t dEdxMIPLog = TMath::Log(2.62463815285237434); //TMath::Log(2.39761562607903311); // [keV/cm] for dX = 2 cm + // Double_t dEdxMIPLog = TMath::Log(2.62463815285237434) + 0.0174824 + 3.57110e-03; ; //TMath::Log(2.39761562607903311); // [keV/cm] for dX = 2 cm + Double_t dEdxMIPLog = TMath::Log(2.62463815285237434); ; //TMath::Log(2.39761562607903311); // [keV/cm] for dX = 2 cm Double_t MIPBetaGamma10 = TMath::Log10(4.); - // log2dx, charge - Double_t pars[3] = { 1.0, 1.0}; + // log2dx, charge mass + Double_t pars[3] = { 1.0, 1.0, 0.13956995}; Double_t dEdxLog = zMP(&MIPBetaGamma10, pars); fLogkeVperElectron = dEdxMIPLog - dEdxLog; cout << "StdEdxModel:: set scale = " << Form("%5.1f",1e3*keVperElectron()) << " eV/electron" << endl; @@ -78,7 +95,16 @@ StdEdxModel::StdEdxModel() : mdNdx(0), fScale(1) StdEdxModel::~StdEdxModel() { fgStdEdxModel = 0; SafeDelete(mdNdx); + SafeDelete(mdNdxL10); + SafeDelete(mLndNdxL10); + SafeDelete(mLndNdxL10Smooth); + SafeDelete(mLndNdxL10Spline5); SafeDelete(fGGaus); + SafeDelete(fGausExp); + SafeDelete(fpol2F); + SafeDelete(fpol5F); + SafeDelete(fpol6F); + SafeDelete(fpol7F); } //________________________________________________________________________________ Double_t StdEdxModel::dNdx(Double_t poverm, Double_t charge) { @@ -93,8 +119,40 @@ Double_t StdEdxModel::dNdx(Double_t poverm, Double_t charge) { Double_t w3 = 121.4139*w2 + 0.0378*TMath::Sin(190.7165*w2); Q_eff *= 1. -w1*TMath::Exp(-w3); } - if (mdNdx) return fScale*Q_eff*Q_eff*mdNdx->Interpolate(poverm); - return 0; + Double_t dNdx = 0; + if ( Stspline3LndNdxL10::instance() && Stspline3LndNdxL10::instance()->IsValid()) { + dNdx = TMath::Exp(Stspline3LndNdxL10::instance()->Func()->Eval(TMath::Log10(poverm))); + } else {// old + if (mLndNdxL10Spline5) { + dNdx = TMath::Exp(mLndNdxL10Spline5->Eval(TMath::Log10(poverm))); + } else if (mLndNdxL10Smooth) { + dNdx = TMath::Exp(mLndNdxL10Smooth->Interpolate(TMath::Log10(poverm))); + } else if (mLndNdxL10) { + dNdx = TMath::Exp(mLndNdxL10->Interpolate(TMath::Log10(poverm))); + } else if (mdNdxL10) { + dNdx = mdNdxL10->Interpolate(TMath::Log10(poverm)); + } else if (mdNdx) { + dNdx = mdNdx->Interpolate(poverm); + } + } + return fScale*Q_eff*Q_eff*dNdx; +} +//________________________________________________________________________________ +Double_t StdEdxModel::dNdxL10func(Double_t *x, Double_t *p) { + Double_t bgL10 = x[0]; + Double_t bg = TMath::Power(10., bgL10); + return instance()->dNdx(bg, 1.); +} +//________________________________________________________________________________ +TF1 *StdEdxModel::dNdxL10F() { + TF1 *f = new TF1("dNdxL10F",dNdxL10func,-2,5,0); + return f; +} +//________________________________________________________________________________ +Double_t StdEdxModel::dNdE() { + static Double_t cLog10 = TMath::Log(10.); + Double_t dE = TMath::Exp(cLog10*mdNdEL10->GetRandom()); + return dE; } //________________________________________________________________________________ Double_t StdEdxModel::gausw(Double_t *x, Double_t *p) { @@ -189,7 +247,7 @@ Double_t StdEdxModel::ggaus(Double_t *x, Double_t *p) { ksi = mu - w*m_0; // Double_t mean = ksi + w * muz; } - Double_t par[4] = {NormL, ksi, w, alpha}; + Double_t par[5] = {NormL, ksi, w, alpha, 0}; Double_t V = gausw(x, par); if (der) { #if 0 /* Derivatives */ @@ -398,7 +456,7 @@ TF1 *StdEdxModel::FParam(Int_t l) { } //________________________________________________________________________________ Double_t StdEdxModel::Prob(Double_t /* log(nE/Np) */ ee, Double_t Np, Double_t *der) { // GG: ggaus - Double_t params[4] = {0}; + Double_t params[5] = {0}; Double_t V = 0; if (! der) { Parameters(Np, ¶ms[1]); @@ -473,7 +531,7 @@ TF1 *StdEdxModel::FProbDer() { #if 0 //________________________________________________________________________________ Double_t StdEdxModel::ProbEx(Double_t /* log(nE/Np) */ ee, Double_t Np, Double_t *der) { // GEX : gausexp - Double_t params[4] = {0}; + Double_t params[5] = {0}; Double_t V = 0; if (! der) { Parameters(Np, ¶ms[1]); @@ -496,12 +554,14 @@ Double_t StdEdxModel::ProbdEGeVlog(Double_t dEGeVLog, Double_t Np, Double_t *der } //________________________________________________________________________________ Double_t StdEdxModel::zMP(Double_t *x, Double_t *p) { // log(keV/cm) - Double_t log10bg = x[0]; - Double_t pOverMRC = TMath::Power(10., log10bg); + Double_t bgL10 = x[0]; + Double_t pOverMRC = TMath::Power(10., bgL10); + // Double_t bgL10 = TMath::Log10(pOverM); Double_t log2dx = p[0]; Double_t charge = p[1]; + Double_t mass = p[2]; Double_t dx = TMath::Power( 2., log2dx); - Double_t dNdx = StdEdxModel::instance()->dNdxEff(pOverMRC, charge); // */dNdxVsBgC*.root [-1.5,5] + Double_t dNdx = StdEdxModel::instance()->dNdxEff(pOverMRC, charge, mass); // */dNdxVsBgC*.root [-1.5,5] Double_t Np = dNdx*dx; // Double_t NpLog = TMath::Log(Np); // Double_t mu = instance()->Parameter(Np, 0); @@ -510,19 +570,70 @@ Double_t StdEdxModel::zMP(Double_t *x, Double_t *p) { // log(keV/cm) // Double_t dEkeVLog = NpLog + mu -3.13746587897608142e+00 +1.78334647296254700e-01;// + 7.02725079814016507e+00;// - 3.13746587897608142e+00;// 43.4 eV/conducting electron Double_t dEkeVLog = instance()->MukeV(Np); // Parameter(Np, 0); Double_t dEdxLog = dEkeVLog - TMath::Log(dx); - return dEdxLog; + Double_t dEdxCor = 0; +#if 0 + StElectonsDEV_dEdx *EL = StElectonsDEV_dEdx::instance(); + if ( EL && EL->IsValid() && EL->InRange(bgL10)) { + dEdxCor = StElectonsDEV_dEdx::instance()->Func()->Eval(bgL10); + static TF1 *elCor1 = 0; + if (! elCor1) { + Double_t pars1[4] = {-0.2100929, 0.1500455, -0.02743834, 0.001894849}; //electrons [2.1,3.3] + Double_t pars2[4] = {-0.04352509, 0.03437069, -0.002851349, -0.0003568138}; //electornsD + Double_t pars[4] = {0}; + for (Int_t i = 0; i < 4; i++) pars[i] = pars1[i] + pars2[i]; + elCor1 = new TF1("dEdxElCor1","pol3",2.0,3.5); elCor1->SetParameters(pars); + } + dEdxCor += elCor1->Eval(bgL10); + dEdxCor += -7.63891e-03 - 3.57110e-03 ; + } else { // pions and protons + static Double_t pionM = 0.13956995; + static Double_t protonM = 0.9382723; + static Double_t mPionL10 = TMath::Log10(pionM); + static Double_t mProtonL10 = TMath::Log10(protonM); + static Double_t dML10 = mProtonL10 - mPionL10; + Double_t dEdxCorPion = 0; + Double_t dEdxCorProton = 0; + StPionDEV_dEdx *PI = StPionDEV_dEdx::instance(); + if ( PI && PI->IsValid() && PI->InRange(bgL10)) { + dEdxCorPion = PI->Func()->Eval(bgL10); + static TF1 *piCor1 = 0; + if (! piCor1) { + Double_t pars[4] = {0.008567411, 0.01817216, -0.004932309, -0.001434306}; //pionD + piCor1 = new TF1("dEdxPiCor1","pol3",-0.2,1.6); piCor1->SetParameters(pars); + } + dEdxCorPion += piCor1->Eval(bgL10); + } + StProtonDEV_dEdx *P = StProtonDEV_dEdx::instance(); + if (P &&P->IsValid() &&P->InRange(bgL10)) { + dEdxCorProton = P->Func()->Eval(bgL10); + static TF1 *protonCor1 = 0; + if (! protonCor1) { + Double_t pars[6] = {0.01745018 + 0.005654033 - 3.57110e-03 , 0.005726225 -0.00228347, 0.004416636, -0.02814983, 0.1824491, -0.2114645}; //protonD + // Double_t pars2[2] = {0.005654033, -0.00228347}; //proton dEdxG + protonCor1 = new TF1("dEdxProtonCor1","pol5",-0.5,0.8); protonCor1->SetParameters(pars); + } + dEdxCorProton += protonCor1->Eval(bgL10); + } + Double_t mL10 = TMath::Log10(mass); + dEdxCor += dEdxCorPion + (dEdxCorProton - dEdxCorPion)*(mL10 - mPionL10)/dML10; + } + // dEdxCor = 1.66944e-02; + // dEdxCor += 0.0174824; +#endif + return dEdxLog + dEdxCor; } //________________________________________________________________________________ -TF1 *StdEdxModel::ZMP(Double_t log2dx) { +TF1 *StdEdxModel::ZMP(Double_t log2dx, Double_t charge, Double_t mass) { TString fName(Form("New%i",(int)(2*(log2dx+2)))); TF1 *f = (TF1 *) gROOT->GetListOfFunctions()->FindObject(fName); if (! f) { - f = new TF1(fName,zMP,-2,5,2); - f->SetParName(0,"log2dx"); + f = new TF1(fName,zMP,-2,5,3); f->SetLineStyle(4); + f->SetParNames("log2dx","charge","mass"); f->SetParameter(0,log2dx); - f->SetParameter(1, 1.0); // charge - cout << "Create ZMPNew with name " << f->GetName() << " for log2dx = " << log2dx << endl; + f->SetParameter(1, charge); + f->SetParameter(2, mass); + cout << "Create ZMPNew with name " << f->GetName() << " for log2dx = " << log2dx << " and mass = " << mass << endl; } return f; } @@ -628,130 +739,42 @@ TH1D *StdEdxModel::protonEff() { //========= Macro generated from object: Func/ //========= by ROOT version5.34/39 - TH1D *Func__1 = new TH1D("protonEff","",100,-2.0,0.0); - Func__1->SetBinContent(1,0.829974); - Func__1->SetBinContent(2,0.833488); - Func__1->SetBinContent(3,0.836679); - Func__1->SetBinContent(4,0.839568); - Func__1->SetBinContent(5,0.842227); - Func__1->SetBinContent(6,0.844577); - Func__1->SetBinContent(7,0.846659); - Func__1->SetBinContent(8,0.848506); - Func__1->SetBinContent(9,0.850145); - Func__1->SetBinContent(10,0.851604); - Func__1->SetBinContent(11,0.857514); - Func__1->SetBinContent(12,0.863548); - Func__1->SetBinContent(13,0.869116); - Func__1->SetBinContent(14,0.874097); - Func__1->SetBinContent(15,0.878448); - Func__1->SetBinContent(16,0.88217); - Func__1->SetBinContent(17,0.885305); - Func__1->SetBinContent(18,0.887923); - Func__1->SetBinContent(19,0.890128); - Func__1->SetBinContent(20,0.891974); - Func__1->SetBinContent(21,0.893513); - Func__1->SetBinContent(22,0.894796); - Func__1->SetBinContent(23,0.895866); - Func__1->SetBinContent(24,0.896761); - Func__1->SetBinContent(25,0.897513); - Func__1->SetBinContent(26,0.898147); - Func__1->SetBinContent(27,0.898683); - Func__1->SetBinContent(28,0.899139); - Func__1->SetBinContent(29,0.899531); - Func__1->SetBinContent(30,0.899868); - Func__1->SetBinContent(31,0.90016); - Func__1->SetBinContent(32,0.900415); - Func__1->SetBinContent(33,0.900639); - Func__1->SetBinContent(34,0.900836); - Func__1->SetBinContent(35,0.90101); - Func__1->SetBinContent(36,0.901165); - Func__1->SetBinContent(37,0.901304); - Func__1->SetBinContent(38,0.901431); - Func__1->SetBinContent(39,0.90164); - Func__1->SetBinContent(40,0.901817); - Func__1->SetBinContent(41,0.901969); - Func__1->SetBinContent(42,0.902099); - Func__1->SetBinContent(43,0.902212); - Func__1->SetBinContent(44,0.902309); - Func__1->SetBinContent(45,0.902392); - Func__1->SetBinContent(46,0.902469); - Func__1->SetBinContent(47,0.902536); - Func__1->SetBinContent(48,0.902595); - Func__1->SetBinContent(49,0.902646); - Func__1->SetBinContent(50,0.902692); - Func__1->SetBinContent(51,0.902733); - Func__1->SetBinContent(52,0.902769); - Func__1->SetBinContent(53,0.902801); - Func__1->SetBinContent(54,0.90283); - Func__1->SetBinContent(55,0.902857); - Func__1->SetBinContent(56,0.902881); - Func__1->SetBinContent(57,0.902902); - Func__1->SetBinContent(58,0.902922); - Func__1->SetBinContent(59,0.90294); - Func__1->SetBinContent(60,0.902956); - Func__1->SetBinContent(61,0.90297); - Func__1->SetBinContent(62,0.902984); - Func__1->SetBinContent(63,0.902996); - Func__1->SetBinContent(64,0.903007); - Func__1->SetBinContent(65,0.903017); - Func__1->SetBinContent(66,0.903026); - Func__1->SetBinContent(67,0.903034); - Func__1->SetBinContent(68,0.903042); - Func__1->SetBinContent(69,0.903049); - Func__1->SetBinContent(70,0.903055); - Func__1->SetBinContent(71,0.903061); - Func__1->SetBinContent(72,0.903066); - Func__1->SetBinContent(73,0.903071); - Func__1->SetBinContent(74,0.903075); - Func__1->SetBinContent(75,0.903079); - Func__1->SetBinContent(76,0.90308); - Func__1->SetBinContent(77,0.90308); - Func__1->SetBinContent(78,0.90308); - Func__1->SetBinContent(79,0.90308); - Func__1->SetBinContent(80,0.90308); - Func__1->SetBinContent(81,0.90308); - Func__1->SetBinContent(82,0.90308); - Func__1->SetBinContent(83,0.90308); - Func__1->SetBinContent(84,0.90308); - Func__1->SetBinContent(85,0.90308); - Func__1->SetBinContent(86,0.90308); - Func__1->SetBinContent(87,0.90308); - Func__1->SetBinContent(88,0.90308); - Func__1->SetBinContent(89,0.90308); - Func__1->SetBinContent(90,0.90308); - Func__1->SetBinContent(91,0.90308); - Func__1->SetBinContent(92,0.90308); - Func__1->SetBinContent(93,0.90308); - Func__1->SetBinContent(94,0.90308); - Func__1->SetBinContent(95,0.90308); - Func__1->SetBinContent(96,0.90308); - Func__1->SetBinContent(97,0.90308); - Func__1->SetBinContent(98,0.90308); - Func__1->SetBinContent(99,0.90308); - Func__1->SetBinContent(100,0.90308); - Func__1->SetEntries(700); - Func__1->SetDirectory(0); - Func__1->SetStats(0); - Func__1->SetFillColor(19); - Func__1->SetFillStyle(0); - Func__1->SetLineColor(9); - Func__1->SetLineWidth(3); - Func__1->SetMarkerStyle(20); - Func__1->GetXaxis()->SetTitleOffset(1.2); - // Func__1->Draw(""); - return Func__1; + TH1D *eff = new TH1D("protonEff","",100,-2.0,0.0); + Double_t corr[102] = { 0, + 0.8300, 0.8335, 0.8367, 0.8396, 0.8422, 0.8446, 0.8467, 0.8485, 0.8501, 0.8516, + 0.8575, 0.8635, 0.8691, 0.8741, 0.8784, 0.8822, 0.8853, 0.8879, 0.8901, 0.8920, + 0.8935, 0.8948, 0.8959, 0.8968, 0.8975, 0.8981, 0.8987, 0.8991, 0.8995, 0.8999, + 0.9002, 0.9004, 0.9006, 0.9008, 0.9010, 0.9012, 0.9013, 0.9014, 0.9016, 0.9018, + 0.9020, 0.9021, 0.9022, 0.9023, 0.9024, 0.9025, 0.9025, 0.9026, 0.9026, 0.9027, + 0.9027, 0.9028, 0.9028, 0.9028, 0.9029, 0.9029, 0.9029, 0.9029, 0.9029, 0.9030, + 0.9030, 0.9030, 0.9030, 0.9030, 0.9030, 0.9030, 0.9030, 0.9030, 0.9030, 0.9031, + 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, + 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, + 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, + 0}; + eff->Set(102, corr); + eff->SetDirectory(0); + eff->SetStats(0); + eff->SetFillColor(19); + eff->SetFillStyle(0); + eff->SetLineColor(9); + eff->SetLineWidth(3); + eff->SetMarkerStyle(20); + eff->GetXaxis()->SetTitleOffset(1.2); + // eff->Draw(""); + return eff; } //________________________________________________________________________________ Double_t StdEdxModel::NpCorrection(Double_t betagamma) { Double_t bgL10 = TMath::Log10(betagamma); - bgL10 = TMath::Max(-2.0, TMath::Min(1.0,bgL10)); + bgL10 = TMath::Max(-2.0, TMath::Min(-1e-3,bgL10)); static TH1D *eff = 0; if (! eff) eff = protonEff(); // return 1.03*eff->Interpolate(bgL10); return eff->Interpolate(bgL10); } //________________________________________________________________________________ -Double_t StdEdxModel::dNdxEff(Double_t poverm, Double_t charge) { +Double_t StdEdxModel::dNdxEff(Double_t poverm, Double_t charge, Double_t mass) { if (!fgStdEdxModel) instance(); Double_t bgMC = bgCorrected(poverm); Double_t dNdxMC = dNdx(bgMC, charge); @@ -759,6 +782,16 @@ Double_t StdEdxModel::dNdxEff(Double_t poverm, Double_t charge) { return dNdx; } //________________________________________________________________________________ +Double_t StdEdxModel::dNdxEffL10func(Double_t *x, Double_t *p) { + Double_t bg = TMath::Power(10., x[0]); + return instance()->dNdxEff(bg, 1., 0.13956995); +} +//________________________________________________________________________________ +TF1 *StdEdxModel::dNdxEffL10F() { + TF1 *f = new TF1("dNdxEffL10F",dNdxEffL10func,-2,5,0); + return f; +} +//________________________________________________________________________________ Double_t StdEdxModel::extremevalueG(Double_t *x, Double_t *p) { Double_t normL = p[0]; Double_t mu = p[1]; diff --git a/StRoot/StBichsel/StdEdxModel.h b/StRoot/StBichsel/StdEdxModel.h index 10f8d1169a8..6658cad1028 100644 --- a/StRoot/StBichsel/StdEdxModel.h +++ b/StRoot/StBichsel/StdEdxModel.h @@ -10,6 +10,7 @@ #include "TFile.h" #include "TString.h" #include "TF1.h" +#include "TSpline.h" class StdEdxParamerization { }; class StdEdxModel { @@ -20,14 +21,18 @@ class StdEdxModel { enum EValType {kProb, kdProbdX, kdProbdY}; virtual ~StdEdxModel(); static StdEdxModel* instance(); - TH1D *GetdNdx() {return mdNdx;} // dN/dx versus beta*gamma static Double_t gausw(Double_t *x, Double_t *p); // vesus ksi, w, alpha static Double_t ggaus(Double_t *x, Double_t *p); // versus mu, sigm, alpha static Double_t ggausD(Double_t *x, Double_t *p, Double_t *der = 0); // versus mu, sigm, alpha wth derivatives static Double_t gausexp(Double_t *x, Double_t *p); // versus mu, sigma, k static Double_t gausexpD(Double_t *x, Double_t *p, Double_t *der = 0); // versus mu, sigma, k - Double_t dNdx(Double_t poverm, Double_t charge = 1.0); - Double_t dNdxEff(Double_t poverm, Double_t charge = 1.0); + Double_t dNdx(Double_t poverm, Double_t charge = 1.0); + TF1 *dNdxL10F(); + static Double_t dNdxL10func(Double_t *x, Double_t *p); + Double_t dNdxEff(Double_t poverm, Double_t charge = 1.0, Double_t mass = 0.13956995); + static Double_t dNdxEffL10func(Double_t *x, Double_t *p); + TF1 *dNdxEffL10F(); + Double_t dNdE(); static Double_t saturationFunc(Double_t *x, Double_t *p); // nP saturation versus beta*gamma from TpcRS (nP/dX - dN/dx_model) TF1 *GGaus() {return fGGaus;} TF1 *GausExp() {return fGausExp;} @@ -63,8 +68,8 @@ class StdEdxModel { Double_t ProbdEGeVlog(Double_t dEGeVLog, Double_t Np, Double_t *der = 0); // probability for give log(dE(GeV)) versus Np void SetScale(Double_t scale = 1.0) {fScale = scale;} Double_t dNdxScale() {return fScale;} - static Double_t zMP(Double_t *x, Double_t *p); // most probable log (dE) versus x = log10(p/M) and p[0] = log2dx and p[1] = charge - TF1 *ZMP(Double_t log2dx = 1); + static Double_t zMP(Double_t *x, Double_t *p); // most probable log (dE) versus x = log10(p/M) and p[0] = log2dx, p[1] = charge, and p[2] = mass + TF1 *ZMP(Double_t log2dx = 1, Double_t charge = 1, Double_t mass = 0.1395699); // from 100 keV Tcut (GEXNor.C) void InitPar(); Double_t tmaxL10eV(Double_t betagamma); // eV @@ -75,16 +80,23 @@ class StdEdxModel { StdEdxModel(); static StdEdxModel *fgStdEdxModel; //! last instance static Int_t _debug; - TH1D *mdNdx; // dN/dx versus beta*gamma Double_t fScale; Double_t fTmaxL10eV; - TF1 *fGGaus; - TF1 *fGausExp; - TF1 *fpol2F; - TF1 *fpol5F; - TF1 *fpol6F; - TF1 *fpol7F; + Char_t beg[1]; //! Double_t fLogkeVperElectron; + TH1D *mdNdx; //! + TH1D *mdNdxL10; //! + TH1D *mLndNdxL10; //! + TH1D *mLndNdxL10Smooth; //! + TSpline5 *mLndNdxL10Spline5; //! + TH1D *mdNdEL10; //! + TF1 *fGGaus; //! + TF1 *fGausExp; //! + TF1 *fpol2F; //! + TF1 *fpol5F; //! + TF1 *fpol6F; //! + TF1 *fpol7F; //! + Char_t end[1]; //! ClassDef(StdEdxModel,0) }; #endif diff --git a/StRoot/StBichsel/StdEdxPull.cxx b/StRoot/StBichsel/StdEdxPull.cxx index c048bf19f50..725e8ba6b5d 100644 --- a/StRoot/StBichsel/StdEdxPull.cxx +++ b/StRoot/StBichsel/StdEdxPull.cxx @@ -2,35 +2,35 @@ #include "StdEdxModel.h" #include "StdEdxPull.h" //________________________________________________________________________________ -Double_t StdEdxPull::EvalPred(Double_t betagamma, UChar_t fit, Int_t charge) { - Double_t dedx_expected; - Double_t dx2 = 1; - if (TMath::Abs(charge) > 1) dx2 = TMath::Log2(5.); +Double_t StdEdxPull::EvalPred(Double_t betagamma, UChar_t fit, Int_t charge, Double_t mass) { + Double_t dedx_expected = 0; if (! fit) { // I70 + Double_t dx2 = 1; + if (TMath::Abs(charge) > 1) dx2 = TMath::Log2(5.); dedx_expected = 1.e-6*charge*charge*Bichsel::Instance()->GetI70M(TMath::Log10(betagamma),dx2); } else if ( fit == 1) { // Ifit // dedx_expected = 1.e-6*charge*charge*TMath::Exp(Bichsel::Instance()->GetMostProbableZ(TMath::Log10(betagamma),dx2)); - Double_t par[3] = {TMath::Log10(betagamma), 1., (Double_t) charge}; + Double_t par[4] = {TMath::Log10(betagamma), 1., (Double_t) charge, mass}; dedx_expected = 1.e-6*TMath::Exp(StdEdxModel::instance()->zMP(par,&par[1])); } else { // dNdx_eff // dedx_expected = StdEdxModel::instance()->dNdx(betagamma,charge); - dedx_expected = StdEdxModel::instance()->dNdxEff(betagamma,charge); + dedx_expected = StdEdxModel::instance()->dNdxEff(betagamma,charge, mass); } return dedx_expected; } //________________________________________________________________________________ -Double_t StdEdxPull::EvalDeV(Double_t dEdx, Double_t betagamma, UChar_t fit, Int_t charge) { - return TMath::Log(dEdx/EvalPred(betagamma, fit, charge)); +Double_t StdEdxPull::EvalDeV(Double_t dEdx, Double_t betagamma, UChar_t fit, Int_t charge, Double_t mass) { + return TMath::Log(dEdx/EvalPred(betagamma, fit, charge, mass)); } //________________________________________________________________________________ -Double_t StdEdxPull::Eval(Double_t dEdx, Double_t dEdxError, Double_t betagamma, UChar_t fit, Int_t charge) { - return (dEdxError > 0) ? EvalDeV(dEdx, betagamma, fit, charge)/dEdxError : -999; +Double_t StdEdxPull::Eval(Double_t dEdx, Double_t dEdxError, Double_t betagamma, UChar_t fit, Int_t charge, Double_t mass) { + return (dEdxError > 0) ? EvalDeV(dEdx, betagamma, fit, charge, mass)/dEdxError : -999; } //________________________________________________________________________________ Double_t StdEdxPull::EvalPred2(Double_t betagamma, Double_t dx2, UChar_t fit, Int_t charge) { - Double_t dedx_expected; - if (TMath::Abs(charge) > 1) dx2 = TMath::Log2(5.); + Double_t dedx_expected = 0; if (! fit) { // I70 + if (TMath::Abs(charge) > 1) dx2 = TMath::Log2(5.); dedx_expected = 1.e-6*charge*charge*Bichsel::Instance()->GetI70M(TMath::Log10(betagamma),dx2); } else if ( fit == 1) { // Ifit dedx_expected = 1.e-6*charge*charge*TMath::Exp(Bichsel::Instance()->GetMostProbableZ(TMath::Log10(betagamma),dx2)); diff --git a/StRoot/StBichsel/StdEdxPull.h b/StRoot/StBichsel/StdEdxPull.h index b089aa2cb76..2f85b8ba8be 100644 --- a/StRoot/StBichsel/StdEdxPull.h +++ b/StRoot/StBichsel/StdEdxPull.h @@ -3,9 +3,9 @@ #include "Rtypes.h" #include "TMath.h" namespace StdEdxPull { - Double_t EvalPred(Double_t betagamma, UChar_t fit = 0, Int_t charge=1); - Double_t EvalDeV(Double_t dEdx, Double_t betagamma, UChar_t fit = 0, Int_t charge=1); - Double_t Eval(Double_t dEdx, Double_t dEdxError, Double_t betagamma, UChar_t fit = 0, Int_t charge=1); + Double_t EvalPred(Double_t betagamma, UChar_t fit = 0, Int_t charge=1, Double_t mass = 0); + Double_t EvalDeV(Double_t dEdx, Double_t betagamma, UChar_t fit = 0, Int_t charge=1, Double_t mass = 0); + Double_t Eval(Double_t dEdx, Double_t dEdxError, Double_t betagamma, UChar_t fit = 0, Int_t charge=1, Double_t mass = 0); Double_t EvalPred2(Double_t betagamma, Double_t dx2, UChar_t fit = 0, Int_t charge=1); Double_t EvalDeV2(Double_t dEdx, Double_t betagamma, Double_t dx2, UChar_t fit = 0, Int_t charge=1); Double_t Eval2(Double_t dEdx, Double_t dEdxError, Double_t betagamma, Double_t dx2, UChar_t fit = 0, Int_t charge=1); diff --git a/StRoot/StBichsel/bichselG10.C b/StRoot/StBichsel/bichselG10.C index 08fa1b57370..ec9bb51b4eb 100644 --- a/StRoot/StBichsel/bichselG10.C +++ b/StRoot/StBichsel/bichselG10.C @@ -85,10 +85,27 @@ #include "StBichsel/StdEdxPull.h" #include "TLegend.h" #include "TROOT.h" +#include "TString.h" #else class Bichsel; #endif -Bichsel *m_Bichsel = 0; +// function object (functor) +struct MyDerivFunc { + MyDerivFunc(TF1 * f): fFunc(f) {} + double operator() (double *x, double * ) const { + return fFunc->Derivative(*x); + } + TF1 * fFunc; +}; +//________________________________________________________________________________ +TF1 *Derivative(TF1 *f1) { + MyDerivFunc * deriv = new MyDerivFunc(f1); + TString name("der_"); + name += f1->GetName(); + TF1 * f2 = new TF1(name,deriv, f1->GetXmin(), f1->GetXmax(), 0, "MyDerivFunc"); + return f2; +} +//Bichsel *m_Bichsel = 0; const Int_t NMasses = 20; const Double_t kAu2Gev=0.9314943228; struct Part_t { @@ -146,8 +163,8 @@ Double_t bichselZ(Double_t *x,Double_t *par) { poverm *= charge; dx2 = TMath::Log2(5.); } - return TMath::Log10(scale*charge*charge*TMath::Exp(m_Bichsel->GetMostProbableZ(TMath::Log10(poverm),dx2)));//TMath::Exp(7.81779499999999961e-01)); - // Charge*Charge* (TMath::Exp(m_Bichsel->GetMostProbableZM(TMath::Log10(TMath::Abs(Charge)*p/M),dx2))) + return TMath::Log10(scale*charge*charge*TMath::Exp(Bichsel::Instance()->GetMostProbableZ(TMath::Log10(poverm),dx2)));//TMath::Exp(7.81779499999999961e-01)); + // Charge*Charge* (TMath::Exp(Bichsel::Instance()->GetMostProbableZM(TMath::Log10(TMath::Abs(Charge)*p/M),dx2))) // return TMath::Log10(1e6*scale*StdEdxPull::EvalPred(poverm, 1, charge)); } //________________________________________________________________________________ @@ -158,7 +175,7 @@ Double_t dEdxModelZ(Double_t *x,Double_t *par) { //new dEdxModel if (mass < 0) {mass = - mass; scale = 2;} Double_t poverm = pove/mass; Double_t charge = par[1]; - poverm *= charge; + poverm *= TMath::Abs(charge); return TMath::Log10(1e6*scale*StdEdxPull::EvalPred(poverm, 1, charge)); } //________________________________________________________________________________ @@ -172,7 +189,7 @@ Double_t bichselZM(Double_t *x,Double_t *par) { Double_t dx2 = 1; if (par[1] > 1.0) { charge = par[1]; - poverm *= charge; + poverm *= TMath::Abs(charge); dx2 = TMath::Log2(5.); } return TMath::Log10(1e6*scale*StdEdxPull::EvalPred(poverm, 1, charge)); @@ -224,7 +241,7 @@ Double_t bichsel70Trs(Double_t *x,Double_t *par) { dx2 = TMath::Log2(5.); } scale *= charge*charge; - return TMath::Log10(scale*TMath::Exp(m_Bichsel->I70Trs(part,TMath::Log10(poverm)))); + return TMath::Log10(scale*TMath::Exp(Bichsel::Instance()->I70Trs(part,TMath::Log10(poverm)))); } //________________________________________________________________________________ Double_t bichselZTrs(Double_t *x,Double_t *par) { @@ -243,7 +260,7 @@ Double_t bichselZTrs(Double_t *x,Double_t *par) { dx2 = TMath::Log2(5.); } scale *= charge*charge; - return TMath::Log10(scale*TMath::Exp(m_Bichsel->IfitTrs(part,TMath::Log10(poverm))));//TMath::Exp(7.81779499999999961e-01)); + return TMath::Log10(scale*TMath::Exp(Bichsel::Instance()->IfitTrs(part,TMath::Log10(poverm))));//TMath::Exp(7.81779499999999961e-01)); } #endif //________________________________________________________________________________ @@ -349,17 +366,15 @@ Double_t aleph70(Double_t *x,Double_t *par) { } #endif /* __CINT__ */ //________________________________________________________________________________ -void bichselG10(const Char_t *type="z", Int_t Nhyps = 9, Bool_t rigidity = kFALSE) { +void bichselG10(const Char_t *type="zN", Int_t Nhyps = 9, Bool_t rigidity = kFALSE) { fgRigidity = rigidity; - if (! m_Bichsel) { - if (gClassTable->GetID("StBichsel") < 0) { - gSystem->Load("libTable"); - gSystem->Load("St_base"); - gSystem->Load("StarClassLibrary"); - gSystem->Load("StBichsel"); - } - m_Bichsel = Bichsel::Instance(); + if (gClassTable->GetID("StBichsel") < 0) { + gSystem->Load("libTable"); + gSystem->Load("St_base"); + gSystem->Load("StarClassLibrary"); + gSystem->Load("StBichsel"); } + // m_Bichsel = Bichsel::Instance(); TString Type(type); TLegend *leg = new TLegend(0.85,0.45,0.95,0.9,""); Double_t xmax = 4; diff --git a/StRoot/StBichsel/spline3.h b/StRoot/StBichsel/spline3.h new file mode 100644 index 00000000000..05a55e42c2f --- /dev/null +++ b/StRoot/StBichsel/spline3.h @@ -0,0 +1,33 @@ +/* spline3.h */ +/* This file was made by the idl compiler "stic". Do not edit. +** This was generated for version '(unspecified)' +** Instead, edit the source idl file and then re-run the compiler. +** For help, type contact Craig Tull or Herb Ward. */ +/* COMMENTS FROM IDL FILE: + spine3.idl + + Table: spine3 + + description: keep parameters and knots for TSpline3 creation + + */ +#ifndef SPLINE3_H +#define SPLINE3_H +#define SPLINE3_SPEC \ +"struct spline3 { \ + long nknots; \ + double Xknots[50]; \ + double Yknots[50]; \ + double ValBeg; \ + double ValEnd; \ + octet option[8]; \ +};" +typedef struct spline3_st { + int nknots; /* no. of knots <= 50 */ + double Xknots[50]; /* X of knots */ + double Yknots[50]; /* Y of knots */ + double ValBeg; + double ValEnd; + unsigned char option[8]; /* option for boundary condition */ +} SPLINE3_ST; +#endif /* SPLINE3_H */ diff --git a/StRoot/StDetectorDbMaker/StDetectorDbChairs.cxx b/StRoot/StDetectorDbMaker/StDetectorDbChairs.cxx index 13379580420..efbd2fc4546 100644 --- a/StRoot/StDetectorDbMaker/StDetectorDbChairs.cxx +++ b/StRoot/StDetectorDbMaker/StDetectorDbChairs.cxx @@ -1,6 +1,7 @@ #include #include #include "StDetectorDbMaker.h" +#include "TROOT.h" #include "TEnv.h" #include "TF1.h" #include "TMath.h" @@ -12,49 +13,59 @@ using namespace ROOT::Math; #include "St_db_Maker/St_db_Maker.h" #if 0 #include "tables/St_tpcCorrection_Table.h" +#include "tables/St_pidCorrection_Table.h" #include "tables/St_tpcSectorT0offset_Table.h" #include "tables/St_tofTrayConfig_Table.h" -#define DEBUGTABLED(STRUCT) PrintTable(#STRUCT,table ) -#define makeString(PATH) # PATH -#define CHECKTABLED(C_STRUCT) \ - if (table->InheritsFrom("St_" makeSTRING(C_STRUCT))) { \ +#define makeSTRING(PATH) # PATH +#define CHECKTABLED(C_STRUCT) \ + if (table->InheritsFrom("St_" # C_STRUCT)) { \ St_ ## C_STRUCT *t = (St_ ## C_STRUCT *) table ; \ - ## C_STRUCT ## _st *s = t->GetTable(); Nrows = s->nrows; \ - ## C_STRUCT ## _st def = {0}; \ + C_STRUCT ## _st *s = t->GetTable(); Nrows = s->nrows; \ + C_STRUCT ## _st def = {0}; \ iprt = kFALSE; \ - Int_t shift = 0; \ - Int_t NrowSize = t->GetRowSize(); \ - if (! strcmp(makeSTRING(C_STRUCT),"Survey")) {shift = 4; NrowSize = 12*8;}\ - if (! strcmp(makeSTRING(C_STRUCT),"tpcSectorT0offset")) {for (Int_t i = 0; i < 24; i++) def->t0[i] = -22.257;} \ - if (! strcmp(makeSTRING(C_STRUCT),"tofTrayConfig")) {def->entries = 120; for (Int_t i = 0; i < 120; i++) {def->iTray[i] = i+1; def->nModules[i] = 32;} \ - for (Int_t i = 0; i < table->GetNRows(); i++, s++) { \ - if (memcmp(&def+shift, s+shift, NrowSize)) {iprt = kTRUE; break;} \ - } \ + Int_t shift = 0; \ + Int_t NrowSize = t->GetRowSize(); \ + if (! strcmp(makeSTRING(C_STRUCT),"Survey")) {shift = 4; NrowSize = 12*8;} \ + if (! strcmp(makeSTRING(C_STRUCT),"tpcSectorT0offset")) {for (Int_t i = 0; i < 24; i++) def.t0[i] = -22.257;} \ + if (! strcmp(makeSTRING(C_STRUCT),"tofTrayConfig")) {def.entries = 120; for (Int_t i = 0; i < 120; i++) {def.iTray[i] = i+1; def.nModules[i] = 32;} \ + for (Int_t i = 0; i < table->GetNRows(); i++, s++) { \ + if (memcmp(&def+shift, s+shift, NrowSize)) {iprt = kTRUE; break;} \ + } \ + } +} +#define DEBUGTABLE(STRUCT) \ + TDatime t[2]; \ + Bool_t iprt = kTRUE; \ + Int_t Nrows = -1; \ + if (St_db_Maker::GetValidity(table,t) > 0) { \ + if (table->InheritsFrom("St_tpcCorrection") ) { \ + St_tpcCorrection *tt = (St_tpcCorrection *) table; \ + tpcCorrection_st *s = tt->GetTable(); Nrows = s->nrows; \ + } else if (table->InheritsFrom("St_pidCorrection") ) { \ + St_pidCorrection *tt = (St_pidCorrection *) table; \ + pidCorrection_st *s = tt->GetTable(); Nrows = s->nrows; \ + } \ + LOG_WARN << "St_" << makeSTRING(STRUCT) << "C::instance found table " << table->GetName() \ + << " with NRows = " << Nrows << " in db" \ + << Form("\tValidity:%8i.%06i",t[0].GetDate(),t[0].GetTime()) \ + << Form(" - %8i.%06i",t[1].GetDate(),t[1].GetTime()) << endm; \ + if (Nrows > 10) Nrows = 10; \ + CHECKTABLED(tpcCorrection); \ + CHECKTABLED(tpcHVPlanes); \ + CHECKTABLED(Survey); \ + CHECKTABLED(tpcSectorT0offset); \ + CHECKTABLED(tofTrayConfig); \ + if (iprt) { \ + if (table->GetRowSize() < 512) table->Print(0,Nrows); \ + } else { \ + LOG_WARN << "Default table" << endm; \ + } \ } //___________________Debug Print out _____________________________________________________________ -void PrintTable(const Char_t *str, TTable *table) { - DEBUGTABLE(str); - Bool_t iprt = kTRUE; - if (St_db_Maker::GetValidity(table,t) > 0) { - if (table->InheritsFrom("St_tpcCorrection")) { - St_tpcCorrection *tt = (St_tpcCorrection *) table; - tpcCorrection_st *s = tt->GetTable(); Nrows = s->nrows;} - if (Nrows > 10) Nrows = 10; - CHECKTABLE(tpcCorrection); - CHECKTABLE(tpcHVPlanes); - CHECKTABLE(Survey); - CHECKTABLE(tpcSectorT0offset); - CHECKTABLE(tofTrayConfig); - if (iprt) { - if (table->GetRowSize() < 512) table->Print(0,Nrows); - } else { - LOG_WARN << "Default table" << endm; - } - } -} #endif #include "StarChairDefs.h" static Int_t _debug = 0; + //___________________Calibrations/ftpc_____________________________________________________________ #include "StDetectorDbFTPCGas.h" StDetectorDbFTPCGas* StDetectorDbFTPCGas::fgInstance = 0; @@ -76,7 +87,11 @@ MakeChairInstance(tpcGridLeak,Calibrations/tpc/tpcGridLeak); #include "St_tpcOmegaTauC.h" MakeChairInstance(tpcOmegaTau,Calibrations/tpc/tpcOmegaTau); #include "St_tpcDriftVelocityC.h" +#include "St_starClockOnlC.h" MakeChairInstance(tpcDriftVelocity,Calibrations/tpc/tpcDriftVelocity); +Float_t St_tpcDriftVelocityC::timeBucketPitch(Int_t i) { + return laserDriftVelocityEast(i)/St_starClockOnlC::instance()->Frequency()*1e6; // cm +} #include "St_TpcSecRowCorC.h" ClassImp(St_TpcSecRowCorC); #include "St_TpcSecRowBC.h" @@ -350,6 +365,126 @@ Int_t St_tpcCorrectionC::IsActiveChair() const { } return npar; } +#include "St_pidCorrectionC.h" +MakeChairInstance2(pidCorrection,St_pidCorrectionC,dEdxModel/PiDC/pidCorrection); +//________________________________________________________________________________ +Double_t St_pidCorrectionC::CalcCorrection(Int_t i, Double_t x, Double_t z, Int_t NparMax) const { + pidCorrection_st *cor = ((St_pidCorrection *) Table())->GetTable() + i; + return SumSeries(cor, x, z, NparMax); +} +//________________________________________________________________________________ +Double_t St_pidCorrectionC::SumSeries(pidCorrection_st *cor, Double_t x, Double_t z, Int_t NparMax) const { + Double_t Sum = 0; + if (! cor) return Sum; + Int_t N = TMath::Abs(cor->npar)%100; + if (N == 0) return Sum; + if (NparMax > 0) N = NparMax; + Double_t X = x; + if (X < cor->min) return Sum; + if (X > cor->max) return Sum; + Double_t T0, T1, T2; + switch (cor->type) { + case 1: // Tchebyshev [-1,1] + if (cor->min < cor->max) X = -1 + 2*TMath::Max(0.,TMath::Min(1.,(X - cor->min)/( cor->max - cor->min))); + T0 = 1; + Sum = cor->a[0]*T0; + if (N == 1) break; + T1 = X; + Sum += cor->a[1]*T1; + for (int n = 2; n <= N; n++) { + T2 = 2*X*T1 - T0; + Sum += cor->a[n]*T2; + T0 = T1; + T1 = T2; + } + break; + case 2: // Shifted TChebyshev [0,1] + if (cor->min < cor->max) X = TMath::Max(0.,TMath::Min(1.,(X - cor->min)/( cor->max - cor->min))); + T0 = 1; + Sum = cor->a[0]*T0; + if (N == 1) break; + T1 = 2*X - 1; + Sum += cor->a[1]*T1; + for (int n = 2; n <= N; n++) { + T2 = 2*(2*X - 1)*T1 - T0; + Sum += cor->a[n]*T2; + T0 = T1; + T1 = T2; + } + break; + default: // polynomials + Sum = cor->a[N-1]; + for (int n = N-2; n>=0; n--) Sum = X*Sum + cor->a[n]; + break; + } + return Sum; +} +//________________________________________________________________________________ +Int_t St_pidCorrectionC::IsActiveChair() const { + const St_pidCorrection *tableC = (const St_pidCorrection *) Table(); + Int_t npar = 0; + if (! tableC) return npar; + pidCorrection_st *cor = tableC->GetTable(); + Int_t N = tableC->GetNRows(); + if (! cor || ! N) { + return npar; + } + N = cor->nrows; + for (Int_t i = 0; i < N; i++, cor++) { + if (cor->nrows == 0 && cor->idx == 0) continue; + if (TMath::Abs(cor->npar) > 0 || + TMath::Abs(cor->OffSet) > 1.e-7 || + TMath::Abs(cor->min) > 1.e-7 || + TMath::Abs(cor->max) > 1.e-7) { + npar++; + } + } + return npar; +} +//________________________________________________________________________________ +Double_t St_pidCorrectionC::Correction(Double_t X, Int_t part, Int_t det, Int_t charge, Int_t pull, Int_t varT) { + Int_t N = nrows(); + Double_t Correction = 0; + for (Int_t l = 0; l < N; l++) { + pidCorrection_st *row = Struct(l); + if (part != row->particle) continue; + if (det != row->det) continue; + if (pull != row->pull) continue; + if (varT != row->var) continue; + if (row->charge != 0 && charge != row->charge) continue; + if (row->min < row->max && (X < row->min || X > row->max)) continue; + // Correction += CalcCorrection(l, X); + Correction += SumSeries(row, X); + } + return TMath::Exp(Correction); +} +//________________________________________________________________________________ +Double_t St_pidCorrectionC::func(Double_t *x, Double_t *p) { + Int_t part = p[0]; + Int_t det = p[1]; + Int_t charge = p[2]; + Int_t pull = p[3]; + Int_t varT = p[4]; + return TMath::Log(St_pidCorrectionC::instance()->Correction(x[0],part,det,charge,pull,varT)); +} +//________________________________________________________________________________ +TF1 *St_pidCorrectionC::Func(Int_t part, Int_t det, Int_t charge, Int_t pull, Int_t varT) { + TString fName(Form("pidCor_p%i_d%i_c%i_p%i_v%i",part,det,charge,pull,varT)); + TF1 *f = (TF1 *) gROOT->GetListOfFunctions()->FindObject(fName); + if (! f) { + f = new TF1(fName, St_pidCorrectionC::func, -1.2, 0.9, 5); + if (! varT) + f->GetXaxis()->SetTitle("log_{10} P [GeV/c]"); + else + f->GetXaxis()->SetTitle("log_{10} P/M"); + f->FixParameter(0, part); + f->FixParameter(1, det); + f->FixParameter(2, charge); + f->FixParameter(3, pull); + f->FixParameter(4, varT); + } + return f; +} #include "St_TpcRowQC.h" MakeChairInstance2(tpcCorrection,St_TpcRowQC,Calibrations/tpc/TpcRowQ); #include "St_TpcDriftDistOxygenC.h" @@ -428,13 +563,19 @@ Double_t St_GatingGridC::CalcCorrection(Int_t i, Double_t x) {// drift time in m MakeChairInstance2(tpcCorrection,St_GatingGridBC,Calibrations/tpc/GatingGridB); //________________________________________________________________________________ Double_t St_GatingGridBC::CalcCorrection(Int_t i, Double_t x) {// drift time in microseconds - Double_t value = -10; - if (x <= t0(i)) return value; - Double_t corD = 1. - TMath::Exp(-(x-t0(i))/(settingTime(i)/4.6)); - if (corD < 1e-4) return value; - return TMath::Log(corD); + const tpcCorrection_st *s = Struct(i); + Double_t T0 = s->a[0]; + Double_t SettingTime = s->a[1]; + // Double_t tau = (x-t0(i))/(settingTime(i)/4.6); + Double_t tau = (x - T0)/(SettingTime/4.6); + if (tau > 10) return 0; + if (tau > 0) { + Double_t corD = 1. - TMath::Exp(-tau); + return TMath::Log(corD); + } else { + return -10; + } } - #include "St_TpcCurrentCorrectionC.h" //MakeChairInstance2(tpcCorrection,St_TpcCurrentCorrectionC,Calibrations/tpc/TpcCurrentCorrection); ClassImp(St_TpcCurrentCorrectionC); @@ -2207,14 +2348,15 @@ double StTpcBXT0CorrEPDC::getCorrection (double epdTAC, double driftVelocity, do } #include "St_tpcT0BXC.h" MakeChairInstance(tpcT0BX,Calibrations/tpc/tpcT0BX); -Double_t St_tpcT0BXC::getT0(Double_t values[6]) const { // xxxEarliestTDC (W+E)/2 variables for "vpd","bbc","epd","zdc", "TAC", "CAVdt" +Double_t St_tpcT0BXC::getT0(Double_t values[7]) const { // xxxEarliestTDC (W+E)/2 variables for "vpd","bbc","epd","zdc", "TAC", "CAVdt" Double_t T0 = 0; Double_t T0W = 0, T0WW = 0; static Int_t __debug = 0; if (__debug > 0) { + if (__debug > 1) Table()->Print(0,10); LOG_INFO << "St_tpcT0BXC::getT0:"; } - for (Int_t i = 0; i < 6; i++) { + for (Int_t i = 0; i < 7; i++) { if (values[i] < xmin(i) || values[i] > xmax(i)) continue; if (detId(i) < 0) continue; // if (detId(i) != 6 && (CountPs(i) < 10 || CountPs(i) > 100)) continue; @@ -2241,6 +2383,18 @@ Double_t St_tpcT0BXC::getT0(Double_t values[6]) const { // xxxEarliestTDC (W+E)/ } #include "St_starTriggerDelayC.h" MakeChairOptionalInstance2(starTriggerDelay,St_starTriggerDelayC,Calibrations/tpc/starTriggerDelay); +//________________________________________________________________________________ +Float_t St_starTriggerDelayC::TrigT0(Int_t i) const { + return clocks(i)/(1e-6*St_starClockOnlC::instance()->CurrentFrequency()) + tZero(i); +} +//________________________________________________________________________________ +Float_t St_starTriggerDelayC::TrigT0GG(Int_t io, Int_t i) const { + Float_t delay = TrigT0(i); // usec add Gating Grid cable + delay += 0.150; // length 30 m, delay 150ns for both Inner and Outer sectors, measured by A.Lebedev 10/13/2023 + if (! io) delay += -0.123; + else delay += -0.502; + return delay; +} //__________________Calibrations/trg______________________________________________________________ #include "St_defaultTrgLvlC.h" MakeChairInstance(defaultTrgLvl,Calibrations/trg/defaultTrgLvl); @@ -2352,7 +2506,6 @@ MakeChairInstance2(spaceChargeCor,St_spaceChargeCorR2C,Calibrations/rich/spaceCh #include "St_MagFactorC.h" MakeChairInstance(MagFactor,RunLog/MagFactor); //_________________RunLog/onl_______________________________________________________________ -#include "St_starClockOnlC.h" MakeChairInstance(starClockOnl,RunLog/onl/starClockOnl); //________________________________________________________________________________ starClockOnl_st *St_starClockOnlC::Struct(Int_t i) { @@ -2413,8 +2566,13 @@ Float_t St_beamInfoC::GammaBlue() { return gamma; } //________________________________________________________________________________ -Float_t St_beamInfoC::Gamma() { - return GammaYellow(); +Float_t St_beamInfoC::GammaCMS() { + Double_t eBlue = kProtonMass*GammaBlue(); + Double_t eYellow = kProtonMass*GammaYellow(); + Double_t Etot = eBlue + eYellow; + Double_t Ecm = SqrtS(); + Double_t gammaCM = Etot/Ecm; + return gammaCM; } //________________________________________________________________________________ Float_t St_beamInfoC::BetaBlue() { @@ -2427,8 +2585,12 @@ Float_t St_beamInfoC::BetaYellow() { return -TMath::Sqrt(1 - 1./(gamma*gamma)); } //________________________________________________________________________________ -Float_t St_beamInfoC::Beta() { - return -BetaBlue(); +Float_t St_beamInfoC::BetaCMS() { + Double_t eBlue = kProtonMass*GammaBlue(); + Double_t eYellow = kProtonMass*GammaYellow(); + Double_t Etot = eBlue + eYellow; + Double_t betaCM = (eBlue*BetaBlue() + eYellow*BetaYellow())/Etot; + return betaCM; } //________________________________________________________________________________ Float_t St_beamInfoC::SqrtS() { @@ -2439,12 +2601,8 @@ Float_t St_beamInfoC::SqrtS() { } //________________________________________________________________________________ Float_t St_beamInfoC::Ycms() { - Double_t eBlue = kProtonMass*GammaBlue(); - Double_t eYellow = kProtonMass*GammaYellow(); - Double_t Etot = eBlue + eYellow; - Double_t Ecm = SqrtS(); - Double_t gammaCM = Etot/Ecm; - Double_t betaCM = (eBlue*BetaBlue() + eYellow*BetaYellow())/Etot; + Double_t gammaCM = GammaCMS(); + Double_t betaCM = BetaCMS(); Double_t yCM = TMath::Log(gammaCM*(1 + betaCM)); return yCM; } @@ -2452,7 +2610,7 @@ Float_t St_beamInfoC::Ycms() { Float_t St_beamInfoC::Frequency() { static Double_t l = 3833.845*9.99999253982746361e-01;// *9.99988614393081399e-01;// *9.99998896969437556e-01; // RHIC perimetr static Double_t NB = 120; // no. of buches, can be changed - Double_t frequency = 1e-6*NB*Beta()*TMath::C()/l; + Double_t frequency = 1e-6*NB*BetaYellow()*TMath::C()/l; return TMath::Abs(frequency); } //________________________________________________________________________________ @@ -3026,12 +3184,17 @@ St_tofSimResParamsC *St_tofSimResParamsC::instance() { DEBUGTABLE(tofSimResParams); fgInstance = new St_tofSimResParamsC(table); mAverageTimeResTof=0; - for ( int i = 0; i < 120; i++ ){ // nTrays - for ( int j = 0; j < 192; j++ ){ - size_t index = i * 120 + j; - params[i][j] = St_tofSimResParamsC::instance()->resolution()[index]; - mAverageTimeResTof+=params[i][j]; - LOG_DEBUG << "tray:" << i << ", mod cell:" << j << " = " << St_tofSimResParamsC::instance()->resolution()[index] << " == " << params[i][j] << endm; + const tofSimResParams_st *row = fgInstance->Struct(0); + for ( Int_t i = 0; i < 120; i++ ){ // nTrays + for ( Int_t j = 0; j < 192; j++ ){ + Int_t index = i * 120 + j; + params[i][j] = row->resolution[index]; + mAverageTimeResTof += params[i][j]; +#if 0 + if (_debug) { + LOG_DEBUG << "tray:" << i << ", mod cell:" << j << " = " << row->resolution[index] << " == " << params[i][j] << endm; + } +#endif } } mAverageTimeResTof=mAverageTimeResTof/(120*192); @@ -3050,11 +3213,8 @@ Double_t St_tofSimResParamsC::timeres_tof(UInt_t itray, UInt_t imodule, UInt_t i * then returns a single vertex resolution (in ps) for use in embedding w/ vpdStart */ Double_t result = 8.5e-11; - if ( itray > 120 || imodule > 32 || icell > 6 ) - return result; - + if ( itray > 120 || imodule > 32 || icell > 6 ) return result; return params[ itray ][ imodule * 6 + icell ]; - } #include "St_vpdDelayC.h" MakeChairInstance(vpdDelay,Calibrations/tof/vpdDelay); diff --git a/StRoot/StDetectorDbMaker/StDetectorDbMaker.cxx b/StRoot/StDetectorDbMaker/StDetectorDbMaker.cxx old mode 100755 new mode 100644 diff --git a/StRoot/StDetectorDbMaker/St_beamInfoC.h b/StRoot/StDetectorDbMaker/St_beamInfoC.h index 63a2f5ed519..8b7e60d0b78 100644 --- a/StRoot/StDetectorDbMaker/St_beamInfoC.h +++ b/StRoot/StDetectorDbMaker/St_beamInfoC.h @@ -42,12 +42,12 @@ class St_beamInfoC : public TChair { Float_t getYellowBunchIntensity(Int_t i=0){return yellowBunchIntensity(i);} Float_t getYellowFillNumber(Int_t i=0) {return yellowFillNumber(i);} Bool_t IsFixedTarget(); - Float_t Gamma(); + Float_t GammaCMS(); Float_t GammaYellow(); Float_t GammaBlue(); Float_t BetaYellow(); Float_t BetaBlue(); - Float_t Beta(); + Float_t BetaCMS(); Float_t Frequency(); Float_t SqrtS(); Float_t Ycms(); diff --git a/StRoot/StDetectorDbMaker/St_pidCorrectionC.h b/StRoot/StDetectorDbMaker/St_pidCorrectionC.h new file mode 100644 index 00000000000..56f16b8de6e --- /dev/null +++ b/StRoot/StDetectorDbMaker/St_pidCorrectionC.h @@ -0,0 +1,55 @@ +#ifndef St_pidCorrectionC_h +#define St_pidCorrectionC_h + +#include "TChair.h" +#include "tables/St_pidCorrection_Table.h" +#include "StEvent/StPidParticleDefinition.h" +#include "StEvent/StEnumerations.h" +#include "TF1.h" +class St_pidCorrectionC : public TChair { + public: + enum PiDStatusIDs { // from StdEdxY2Maker/StTrackCombPiD.h + kUndef = kUndefinedMethodId, + kI70 = kTruncatedMeanId, + kI70U = kEnsembleTruncatedMeanId, + kFit = kLikelihoodFitId, + kFitU = kWeightedTruncatedMeanId, + kdNdx = kOtherMethodId, + kdNdxU = kOtherMethodId2, + kBTof, kETof, kMtd, kBEmc, kTotal + }; + static St_pidCorrectionC* instance(); + pidCorrection_st *Struct(Int_t i = 0) const {return ((St_pidCorrection*) Table())->GetTable()+i;} + UInt_t getNumRows() const {return GetNRows();} + Int_t idx(Int_t i = 0) const {return Struct(i)->idx;} // row index > 0 if it is real + Int_t nrows(Int_t i = 0) const {return Struct(i)->nrows;} // total no. of real rows in the table; For Db interface (where nrows = 50) + Int_t type(Int_t i = 0) const {return Struct(i)->type;} // type = 0 polymonical fit, use only [min,max] + Int_t var(Int_t i = 0) const {return Struct(i)->var;} // fit variable: 0 => pmomL10, 1 => bgL10 + Int_t particle(Int_t i = 0) const {return Struct(i)->particle;} // from StEvent/StPidParticleDefinition.h : kUndef = -1, kPidElectron = 0, Proton = 1, + // Kaon = 2, Pion = 3, Muon = 4, Deuteron = 5, Triton = 6, + // He3 = 7, Alpha = 8, He6 = 9, Li5 = 10, Li6,= 11, Li7 = 12, Be7 = 13, Be9 = 14, Be10 = 15, B11 = 16 + Int_t charge(Int_t i = 0) const {return Struct(i)->charge;} // +/-1, 0 does not matter + Int_t pull(Int_t i = 0) const {return Struct(i)->pull;} // != 0 calculated pull correction, == 0 to value + Int_t det(Int_t i = 0) const {return Struct(i)->det;} // from StdEdxY2Maker/StTrackPiD.h : kUndef = 0, kI70 = 1, kI70U = 2, kFit = 3, kFitU = 4, kdNdx = 5, + // kdNdxU = 6, kBTof -7 , kETof = 8, kMtd = 9, kBEmc = 10 + Int_t npar(Int_t i = 0) const {return Struct(i)->npar;} // npar < 0, X = exp(x) paramterization; abs(npar) >= 100 cut on range [min.max] + Double_t OffSet(Int_t i = 0) const {return Struct(i)->OffSet;} // for future use + Double_t min(Int_t i = 0) const {return Struct(i)->min;} // fit range + Double_t max(Int_t i = 0) const {return Struct(i)->max;} // + Double_t* a(Int_t i = 0) const {return Struct(i)->a;} // a[npar] + Char_t* comment(Int_t i = 0) const {return Struct(i)->comment;} + Double_t CalcCorrection(Int_t i, Double_t x, Double_t z = 0, Int_t NparMax = -1) const; + Double_t SumSeries(pidCorrection_st *cor, Double_t x, Double_t z = 0, Int_t NparMax = -1) const; + Double_t Correction(Double_t X, Int_t part = kPidPion, Int_t det = kFit, Int_t charge = 0, Int_t pull = 0, Int_t varT = 0); + static Double_t func(Double_t *x, Double_t *p); + TF1* Func(Int_t part = kPidPion, Int_t det = kFit, Int_t charge = 0, Int_t pull = 0, Int_t varT = 0); + Int_t IsActiveChair() const; + protected: + St_pidCorrectionC(St_pidCorrection *table=0) : TChair(table) {} + virtual ~St_pidCorrectionC() {fgInstance = 0;} + private: + static St_pidCorrectionC* fgInstance; + ClassDefChair(St_pidCorrection, pidCorrection_st ) + ClassDef(St_pidCorrectionC,1) //C++ TChair for pidCorrection table class +}; +#endif diff --git a/StRoot/StDetectorDbMaker/St_starTriggerDelayC.h b/StRoot/StDetectorDbMaker/St_starTriggerDelayC.h index 843647a6d00..aab06bac1e0 100644 --- a/StRoot/StDetectorDbMaker/St_starTriggerDelayC.h +++ b/StRoot/StDetectorDbMaker/St_starTriggerDelayC.h @@ -11,6 +11,8 @@ class St_starTriggerDelayC : public TChair { UInt_t getNumRows() const {return GetNRows();} Float_t clocks(Int_t i = 0) const {return Struct(i)->clocks;} Float_t tZero(Int_t i = 0) const {return Struct(i)->tZero;} + Float_t TrigT0(Int_t i = 0) const; // usec + Float_t TrigT0GG(Int_t io = 0, Int_t i = 0) const; // usec add cables protected: St_starTriggerDelayC(St_starTriggerDelay *table=0) : TChair(table) {} virtual ~St_starTriggerDelayC() {if (Table()->IsMarked()) delete GetThisTable(); fgInstance = 0;} diff --git a/StRoot/StDetectorDbMaker/St_tpcDriftVelocityC.h b/StRoot/StDetectorDbMaker/St_tpcDriftVelocityC.h index baa2be3cd6f..e50241a0194 100644 --- a/StRoot/StDetectorDbMaker/St_tpcDriftVelocityC.h +++ b/StRoot/StDetectorDbMaker/St_tpcDriftVelocityC.h @@ -13,6 +13,7 @@ class St_tpcDriftVelocityC : public TChair { Float_t laserDriftVelocityWest(Int_t i = 0) {return Struct(i)->laserDriftVelocityWest;} Float_t cathodeDriftVelocityEast(Int_t i = 0) {return Struct(i)->cathodeDriftVelocityEast;} Float_t cathodeDriftVelocityWest(Int_t i = 0) {return Struct(i)->cathodeDriftVelocityWest;} + Float_t timeBucketPitch(Int_t i = 0); // cm #if 0 Float_t scaleY(Int_t i = 0) {return Struct(i)->scaleY;} #endif diff --git a/StRoot/StDetectorDbMaker/St_tpcPadConfigC.h b/StRoot/StDetectorDbMaker/St_tpcPadConfigC.h index b55e690f0b4..21ee90bcc4e 100644 --- a/StRoot/StDetectorDbMaker/St_tpcPadConfigC.h +++ b/StRoot/StDetectorDbMaker/St_tpcPadConfigC.h @@ -10,58 +10,58 @@ class St_tpcPadConfigC : public TChair { tpcPadConfig_st *Struct(Int_t i=0); UInt_t getNumRows(); UChar_t *itpc(Int_t i=0) {return (((St_tpcPadConfig*) Table())->GetTable(i))->itpc;} - UChar_t iTpc(Int_t sector); - UChar_t iTPC(Int_t sector) {return iTpc(sector);} - Int_t padRows(Int_t sector); - Int_t innerPadRows(Int_t sector); - Int_t innerPadRows48(Int_t sector); - Int_t innerPadRows52(Int_t sector); - Int_t outerPadRows(Int_t sector); - Int_t superInnerPadRows(Int_t sector); - Int_t superOuterPadRows(Int_t sector); - Double_t innerSectorPadWidth(Int_t sector); - Double_t innerSectorPadLength(Int_t sector); - Double_t innerSectorPadPitch(Int_t sector); - Double_t innerSectorRowPitch1(Int_t sector); - Double_t innerSectorRowPitch2(Int_t sector); - Double_t firstPadRow(Int_t sector); - Double_t firstOuterSectorPadRow(Int_t sector); - Double_t lastOuterSectorPadRow(Int_t sector); - Double_t firstRowWidth(Int_t sector); - Double_t lastRowWidth(Int_t sector); - Double_t outerSectorPadWidth(Int_t sector); - Double_t outerSectorPadLength(Int_t sector); - Double_t outerSectorPadPitch(Int_t sector); - Double_t outerSectorRowPitch(Int_t sector); - Double_t outerSectorLength(Int_t sector); - Double_t ioSectorSeparation(Int_t sector); - Double_t innerSectorEdge(Int_t sector); - Double_t outerSectorEdge(Int_t sector); - Double_t innerSectorPadPlaneZ(Int_t sector); - Double_t outerSectorPadPlaneZ(Int_t sector); - Int_t* innerPadsPerRow(Int_t sector); - Int_t* outerPadsPerRow(Int_t sector); - Int_t padsPerRow(Int_t sector, Int_t row = 1); - Double_t* innerRowRadii(Int_t sector); - Double_t* outerRowRadii(Int_t sector); + UChar_t iTpc(Int_t sector = 1); + UChar_t iTPC(Int_t sector = 1) {return iTpc(sector);} + Int_t padRows(Int_t sector = 1); + Int_t innerPadRows(Int_t sector = 1); + Int_t innerPadRows48(Int_t sector = 1); + Int_t innerPadRows52(Int_t sector = 1); + Int_t outerPadRows(Int_t sector = 1); + Int_t superInnerPadRows(Int_t sector = 1); + Int_t superOuterPadRows(Int_t sector = 1); + Double_t innerSectorPadWidth(Int_t sector = 1); + Double_t innerSectorPadLength(Int_t sector = 1); + Double_t innerSectorPadPitch(Int_t sector = 1); + Double_t innerSectorRowPitch1(Int_t sector = 1); + Double_t innerSectorRowPitch2(Int_t sector = 1); + Double_t firstPadRow(Int_t sector = 1); + Double_t firstOuterSectorPadRow(Int_t sector = 1); + Double_t lastOuterSectorPadRow(Int_t sector = 1); + Double_t firstRowWidth(Int_t sector = 1); + Double_t lastRowWidth(Int_t sector = 1); + Double_t outerSectorPadWidth(Int_t sector = 1); + Double_t outerSectorPadLength(Int_t sector = 1); + Double_t outerSectorPadPitch(Int_t sector = 1); + Double_t outerSectorRowPitch(Int_t sector = 1); + Double_t outerSectorLength(Int_t sector = 1); + Double_t ioSectorSeparation(Int_t sector = 1); + Double_t innerSectorEdge(Int_t sector = 1); + Double_t outerSectorEdge(Int_t sector = 1); + Double_t innerSectorPadPlaneZ(Int_t sector = 1); + Double_t outerSectorPadPlaneZ(Int_t sector = 1); + Int_t* innerPadsPerRow(Int_t sector = 1); + Int_t* outerPadsPerRow(Int_t sector = 1); + Int_t padsPerRow(Int_t sector = 1, Int_t row = 1); + Double_t* innerRowRadii(Int_t sector = 1); + Double_t* outerRowRadii(Int_t sector = 1); // taken from StRItpcPadPlane - Int_t numberOfRows(Int_t sector); - Int_t numberOfInnerRows(Int_t sector); - Int_t numberOfInnerRows48(Int_t sector); - Int_t numberOfInnerRows52(Int_t sector); - Int_t numberOfOuterRows(Int_t sector); - Bool_t isRowInRange(Int_t sector, Int_t row); - Double_t radialDistanceAtRow(Int_t sector, Int_t row); - Int_t numberOfPadsAtRow(Int_t sector, Int_t row); - Double_t PadWidthAtRow(Int_t sector, Int_t row); - Double_t PadLengthAtRow(Int_t sector, Int_t row); - Double_t PadPitchAtRow(Int_t sector, Int_t row); - Double_t RowPitchAtRow(Int_t sector, Int_t row); - Int_t indexForRowPad(Int_t sector, Int_t row, Int_t pad); - bool isiTpcSector(Int_t sector) { return iTpc(sector) == 1; } - bool isiTpcPadRow(Int_t sector, Int_t row) { return iTpc(sector) && row >= 1 && row <= numberOfInnerRows(sector); } - bool isInnerPadRow(Int_t sector, Int_t row) { return row <= numberOfInnerRows(sector); } - Int_t IsRowInner(Int_t sector, Int_t row) {return (row <= innerPadRows(sector)) ? 1 : 0;} + Int_t numberOfRows(Int_t sector = 1); + Int_t numberOfInnerRows(Int_t sector = 1); + Int_t numberOfInnerRows48(Int_t sector = 1); + Int_t numberOfInnerRows52(Int_t sector = 1); + Int_t numberOfOuterRows(Int_t sector = 1); + Bool_t isRowInRange(Int_t sector = 1, Int_t row = 1); + Double_t radialDistanceAtRow(Int_t sector = 1, Int_t row = 1); + Int_t numberOfPadsAtRow(Int_t sector = 1, Int_t row = 1); + Double_t PadWidthAtRow(Int_t sector = 1, Int_t row = 1); + Double_t PadLengthAtRow(Int_t sector = 1, Int_t row = 1); + Double_t PadPitchAtRow(Int_t sector = 1, Int_t row = 1); + Double_t RowPitchAtRow(Int_t sector = 1, Int_t row = 1); + Int_t indexForRowPad(Int_t sector = 1, Int_t row = 1, Int_t pad = 1); + bool isiTpcSector(Int_t sector = 1) { return iTpc(sector) == 1; } + bool isiTpcPadRow(Int_t sector = 1, Int_t row = 1) { return iTpc(sector) && row >= 1 && row <= numberOfInnerRows(sector); } + bool isInnerPadRow(Int_t sector = 1, Int_t row = 1) { return row <= numberOfInnerRows(sector); } + Int_t IsRowInner(Int_t sector = 1, Int_t row = 1) {return (row <= innerPadRows(sector)) ? 1 : 0;} protected: St_tpcPadConfigC(St_tpcPadConfig *table=0) : TChair(table) {} virtual ~St_tpcPadConfigC() {fgInstance = 0;} diff --git a/StRoot/StDetectorDbMaker/St_tpcT0BXC.h b/StRoot/StDetectorDbMaker/St_tpcT0BXC.h index a506897d679..dfaed2f326c 100644 --- a/StRoot/StDetectorDbMaker/St_tpcT0BXC.h +++ b/StRoot/StDetectorDbMaker/St_tpcT0BXC.h @@ -21,7 +21,7 @@ class St_tpcT0BXC : public TChair { Float_t dslope(Int_t i = 0) const {return Struct(i)->dslope;} Float_t CountPs(Int_t i = 0) const {return Struct(i)->CountPs;} Float_t dCountPs(Int_t i = 0) const {return Struct(i)->dCountPs;} - Double_t getT0(Double_t values[6]) const; + Double_t getT0(Double_t values[7]) const; protected: St_tpcT0BXC(St_tpcT0BX *table=0) : TChair(table) {} virtual ~St_tpcT0BXC() {fgInstance = 0;} diff --git a/StRoot/StDetectorDbMaker/St_trgTimeOffsetC.h b/StRoot/StDetectorDbMaker/St_trgTimeOffsetC.h index 5a14b83e5c1..c0f4c5bc0bf 100644 --- a/StRoot/StDetectorDbMaker/St_trgTimeOffsetC.h +++ b/StRoot/StDetectorDbMaker/St_trgTimeOffsetC.h @@ -15,6 +15,7 @@ class St_trgTimeOffsetC : public TChair { Float_t triggerTimeOffset(Int_t i = 0) {return 1e-6*(mLaser ? laserOffset(i) : offset(i));} // usec Float_t triggerTimeOffsetWest(Int_t i = 0) {return 1e-6*(mLaser ? laserOffsetW(i) : 0);} // usec void SetLaser(Bool_t k = kTRUE) {mLaser = k;} + Int_t IsLaser() { return mLaser;} protected: St_trgTimeOffsetC(St_trgTimeOffset *table=0) : TChair(table), mLaser(kFALSE) {} virtual ~St_trgTimeOffsetC() {fgInstance = 0;} diff --git a/StRoot/StDetectorDbMaker/StiChairs.cxx b/StRoot/StDetectorDbMaker/StiChairs.cxx index 643e48e79ad..15f9503b8ce 100644 --- a/StRoot/StDetectorDbMaker/StiChairs.cxx +++ b/StRoot/StDetectorDbMaker/StiChairs.cxx @@ -2,14 +2,13 @@ #include "StarChairDefs.h" #include "St_db_Maker/St_db_Maker.h" #include "StDetectorDbMaker/StiHitErrorCalculator.h" -//________________________________________________________________________________ -#include "StiHitErrorCalculator.h" ClassImp(StiHitErrorCalculator); //________________________________________________________________________________ void StiHitErrorCalculator::calculateError(Double_t _z, Double_t _eta, Double_t _tanl, Double_t &ecross, Double_t &edip, Double_t fudgeCactor) const { static const Double_t tenMicrons = 1e-3; static const Double_t min2Err = tenMicrons*tenMicrons; static const Double_t max2Err = 1.; + static const Double_t scale = 1.; const Double_t *Coeff = ((StiHitErrorCalculator *) this)->coeff(); #if 0 Double_t dz = (200.-TMath::Abs(_z+100))/100.; // Local z @@ -27,12 +26,12 @@ void StiHitErrorCalculator::calculateError(Double_t _z, Double_t _eta, Double_t Double_t sinCA = TMath::Sin(Phi); if (TMath::Abs(cosCA)<0.01) cosCA=0.01; Double_t tanCA = sinCA/cosCA; - ecross=fudgeCactor*fudgeCactor*(Coeff[0]+Coeff[1]*dz/(cosCA*cosCA) +Coeff[2]*tanCA*tanCA); + ecross=scale*fudgeCactor*fudgeCactor*(Coeff[0]+Coeff[1]*dz/(cosCA*cosCA) +Coeff[2]*tanCA*tanCA); if (ecross< min2Err) ecross = min2Err; if (ecross> max2Err) ecross = max2Err; Double_t tanDip=_tanl; Double_t cosDipInv2=1+tanDip*tanDip; - edip=fudgeCactor*fudgeCactor*(Coeff[3]+Coeff[4]*dz*cosDipInv2+Coeff[5]*tanDip*tanDip); + edip=scale*fudgeCactor*fudgeCactor*(Coeff[3]+Coeff[4]*dz*cosDipInv2+Coeff[5]*tanDip*tanDip); if (edip< min2Err) edip = min2Err; if (edip> max2Err) edip = max2Err; // Temporary hack for Gene. Increase prompt hit errors @@ -51,3 +50,48 @@ MakeChairInstance2(LocalTrackSeedFinder,StiLocalTrackSeedFinderParameters,Calibr MakeChairInstance2(KalmanTrackFitterParameters,StiKalmanTrackFitterParameters,Calibrations/tracker/KalmanTrackFitterParameters); #include "StiKalmanTrackFinderParameters.h" MakeChairInstance2(KalmanTrackFinderParameters,StiKalmanTrackFinderParameters,Calibrations/tracker/KalmanTrackFinderParameters); +#include "StiTpcHitErrorMDF4.h" +//________________________________________________________________________________ +void StiTpcHitErrorMDF4::convert(Double_t _z, Double_t _eta, Double_t _tanl, Double_t AdcL) { + fxx[0] = 1. - TMath::Abs(_z)/207.707; // Z + Double_t y = TMath::Tan(_eta); + fxx[1] = y*y; // tanP**2 + fxx[2] = _tanl*_tanl; // tanL**2 + fxx[3] = AdcL; // AdcL +} +//________________________________________________________________________________ +void StiTpcHitErrorMDF4::calculateError(Double_t _z, Double_t _eta, Double_t _tanl, + Double_t &ecross, Double_t &edip, + Double_t fudgeFactor, Double_t AdcL, + Double_t *dZ, Double_t *dX) { + static const Double_t tenMicrons = 1e-3; + static const Double_t min2Err = tenMicrons*tenMicrons; + static const Double_t max2Err = 1.; + static const Double_t scale = 1.; + convert(_z, _eta, _tanl, AdcL); + Double_t dPadSigmaSQ = Eval( 0, fxx); + Double_t dTimeSigmaSQ = Eval( 2, fxx); + ecross = scale*padPitch() *padPitch() *dPadSigmaSQ * fudgeFactor; + edip = scale*timePitch()*timePitch()*dTimeSigmaSQ * fudgeFactor; + Int_t fail = 0; + if (ecross< min2Err) {ecross = min2Err; fail++;} + if (ecross> max2Err) {ecross = max2Err; fail++;} + if (edip< min2Err) {edip = min2Err; fail++;} + if (edip> max2Err) {edip = max2Err; fail++;} + if (dZ) { + if (fail) *dZ = 0; + else { + Double_t dTime = Eval( 3, fxx); + *dZ = - timePitch()*dTime * TMath::Sign(1., _z); + } + } + if (dX) { + if (fail) *dX = 0; + else { + Double_t dPad = Eval( 1, fxx); + *dX = - padPitch()*dPad; + } + } +} +MakeChairInstance2(MDFCorrection4,StiTpcInnerHitErrorMDF4,Calibrations/tracker/TpcInnerHitErrorMDF4); +MakeChairInstance2(MDFCorrection4,StiTpcOuterHitErrorMDF4,Calibrations/tracker/TpcOuterHitErrorMDF4); diff --git a/StRoot/StDetectorDbMaker/StiTpcHitErrorMDF4.h b/StRoot/StDetectorDbMaker/StiTpcHitErrorMDF4.h new file mode 100644 index 00000000000..5d2ef96dd7d --- /dev/null +++ b/StRoot/StDetectorDbMaker/StiTpcHitErrorMDF4.h @@ -0,0 +1,55 @@ +#ifndef StiTpcHitErrorMDF4_h +#define StiTpcHitErrorMDF4_h + +#include "Sti/StiNodePars.h" +#include "St_MDFCorrection4C.h" +#include "StDetectorDbMaker/St_tpcDriftVelocityC.h" +#include "StDetectorDbMaker/St_tpcPadConfigC.h" + +class StiTpcHitErrorMDF4 : public St_MDFCorrection4C { + public: + virtual void calculateError(Double_t _z, Double_t _eta, Double_t _tanl, + Double_t &ecross, Double_t &edip, + Double_t fudgeFactor = 1, Double_t AdcL = 5.5, Double_t *dZ = 0, Double_t *dX = 0); + virtual void calculateError(const StiNodePars *pars, + Double_t &ecross, Double_t &edip, + Double_t fudgeFactor = 1, Double_t AdcL = 5.5, Double_t *dZ = 0, Double_t *dX = 0) { + calculateError(pars->z(), pars->eta(), pars->tanl(), ecross, edip, fudgeFactor, AdcL, dZ, dX); + } + protected: + StiTpcHitErrorMDF4(St_MDFCorrection4 *table=0) : St_MDFCorrection4C(table) {} + virtual ~StiTpcHitErrorMDF4() {} + private: + void convert(Double_t _z, Double_t _eta, Double_t _tanl, Double_t AdcL); + Double_t fxx[4]; + virtual Double_t padPitch() = 0; + Double_t timePitch() {return St_tpcDriftVelocityC::instance()->timeBucketPitch();} + ClassDefineChair(StiTpcHitErrorMDF4,St_MDFCorrection4, MDFCorrection4_st ) + ClassDef(StiTpcHitErrorMDF4,1) //C++ TChair for MDFCorrection4 table class +}; +//________________________________________________________________________________ +class StiTpcInnerHitErrorMDF4 : public StiTpcHitErrorMDF4 { + public: + static StiTpcInnerHitErrorMDF4 *instance(); + protected: + StiTpcInnerHitErrorMDF4(St_MDFCorrection4 *table=0) : StiTpcHitErrorMDF4(table) {} + virtual ~StiTpcInnerHitErrorMDF4() {fgInstance = 0;} + private: + static StiTpcInnerHitErrorMDF4* fgInstance; + Double_t padPitch() {return St_tpcPadConfigC::instance()->innerSectorPadPitch();} + ClassDef(StiTpcInnerHitErrorMDF4,1) //C++ TChair for MDFCorrection4 table class +}; +//________________________________________________________________________________ +class StiTpcOuterHitErrorMDF4 : public StiTpcHitErrorMDF4 { + public: + static StiTpcOuterHitErrorMDF4 *instance(); + protected: + StiTpcOuterHitErrorMDF4(St_MDFCorrection4 *table=0) : StiTpcHitErrorMDF4(table) {} + virtual ~StiTpcOuterHitErrorMDF4() {fgInstance = 0;} + private: + static StiTpcOuterHitErrorMDF4* fgInstance; + Double_t padPitch() {return St_tpcPadConfigC::instance()->innerSectorPadPitch();} + ClassDef(StiTpcOuterHitErrorMDF4,1) //C++ TChair for MDFCorrection4 table class +}; + +#endif diff --git a/StRoot/StTpcRSMaker/StTpcRSMaker.cxx b/StRoot/StTpcRSMaker/StTpcRSMaker.cxx index 84ebddc4005..291596a9f77 100644 --- a/StRoot/StTpcRSMaker/StTpcRSMaker.cxx +++ b/StRoot/StTpcRSMaker/StTpcRSMaker.cxx @@ -51,6 +51,8 @@ #include "StDetectorDbMaker/St_TpcPadPedRMSC.h" #include "StEventUtilities/StEbyET0.h" #include "StDetectorDbMaker/St_beamInfoC.h" +#include "StDetectorDbMaker/St_GatingGridBC.h" +#include "StDetectorDbMaker/St_starTriggerDelayC.h" #if 0 #include "StParticleTable.hh" #include "StParticleDefinition.hh" @@ -58,6 +60,7 @@ #include "Altro.h" #include "TRVector.h" #include "StBichsel/Bichsel.h" +#include "StBichsel/StdEdxModel.h" #include "StdEdxY2Maker/StTpcdEdxCorrection.h" // g2t tables #include "tables/St_g2t_track_Table.h" @@ -104,14 +107,16 @@ static Double_t fgTriggerT0 = 0; //! TPC trigger T0 (seconds) is sup static Double_t timeBinMin = -0.5; static Double_t timeBinMax = 44.5; //________________________________________________________________________________ -static const Int_t nx[2] = {200,500}; -static const Double_t xmin[2] = {-10., -6}; -static const Double_t xmax[2] = { 10., 44}; +static const Int_t nx[4] = {200,500, 145, 401}; +static const Double_t xmin[4] = {-10., -6, -0.5, -0.5}; +static const Double_t xmax[4] = { 10., 44,144.5, 400.5}; static const Int_t nz = 42; static const Double_t zmin = -210; static const Double_t zmax = -zmin; // io pt -static TProfile2D *hist[5][3] = {0}; +static TProfile2D *hist[4][2] = {0}; +static TProfile2D *PixelHist[2] = {0}; +static TProfile2D *GainHist[2] = {0}; static const Int_t nChecks = 22; static TH1 *checkList[2][22] = {0}; static TString TpcMedium("TPCE_SENSITIVE_GAS"); @@ -146,9 +151,6 @@ StTpcRSMaker::~StTpcRSMaker() { Int_t StTpcRSMaker::Finish() { // SafeDelete(fTree); if (m_SignalSum) {free(m_SignalSum); m_SignalSum = 0;} - SafeDelete(mdNdx); - SafeDelete(mdNdxL10); - SafeDelete(mdNdEL10); for (Int_t io = 0; io < 2; io++) {// Inner/Outer for (Int_t sec = 0; sec < NoOfSectors; sec++) { if (mShaperResponses[io][sec] && !mShaperResponses[io][sec]->TestBit(kNotDeleted)) {SafeDelete(mShaperResponses[io][sec]);} @@ -180,45 +182,28 @@ Int_t StTpcRSMaker::InitRun(Int_t /* runnumber */) { LOG_ERROR << "StTpcRSMaker::InitRun: mCutEle has not been found in GEANT3 for \"" << TpcMedium.Data() << "\" parameters." << "Probably due to missing Set it to default " << mCutEle << endm; } - LOG_INFO << "StTpcRSMaker:: use H.Bichsel model for dE/dx simulation" << endm; - if (! mdNdEL10 || ! mdNdx) { - const Char_t *path = ".:./StarDb/dEdxModel:./StarDb/global/dEdx" - ":./StRoot/StBichsel:$STAR/StarDb/dEdxModel:$STAR/StarDb/global/dEdx:$STAR/StRoot/StBichsel"; - const Char_t *Files[2] = {"dNdE_Bichsel.root","dNdx_Bichsel.root"}; - for (Int_t i = 0; i < 2; i++) { // Inner/Outer - Char_t *file = gSystem->Which(path,Files[i],kReadPermission); - if (! file) Fatal("StTpcRSMaker::Init","File %s has not been found in path %s",Files[i],path); - else Warning("StTpcRSMaker::Init","File %s has been found as %s",Files[i],file); - TFile *pFile = new TFile(file); - if (i == 0) {mdNdEL10 = (TH1D *) pFile->Get("dNdEL10"); assert(mdNdEL10); mdNdEL10->SetDirectory(0);} - if (i == 1) {mdNdx = (TH1D *) pFile->Get("dNdx"); assert(mdNdx); mdNdx->SetDirectory(0);} - delete pFile; - delete [] file; - } - assert(mdNdEL10 && mdNdx); - } // Distortions if (TESTBIT(m_Mode,kdEdxCorr)) { LOG_INFO << "StTpcRSMaker:: use Tpc dE/dx correction from calibaration" << endm; Long_t Mask = -1; // 64 bits CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrection); CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrectionC); + CLRBIT(Mask,StTpcdEdxCorrection::kEdge); CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrectionMDF); CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrection3MDF); CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrection4MDF); CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrection5MDF); CLRBIT(Mask,StTpcdEdxCorrection::kAdcCorrection6MDF); - CLRBIT(Mask,StTpcdEdxCorrection::kEdge); CLRBIT(Mask,StTpcdEdxCorrection::kEtaCorrection); -#if 0 - CLRBIT(Mask,StTpcdEdxCorrection::kdXCorrection); + CLRBIT(Mask,StTpcdEdxCorrection::knTbk); CLRBIT(Mask,StTpcdEdxCorrection::kzCorrectionC); CLRBIT(Mask,StTpcdEdxCorrection::kzCorrection); +#if 0 + CLRBIT(Mask,StTpcdEdxCorrection::kdXCorrection); CLRBIT(Mask,StTpcdEdxCorrection::kTpcPadTBins); CLRBIT(Mask,StTpcdEdxCorrection::kTanL); CLRBIT(Mask,StTpcdEdxCorrection::kAdcI); CLRBIT(Mask,StTpcdEdxCorrection::knPad); - CLRBIT(Mask,StTpcdEdxCorrection::knTbk); CLRBIT(Mask,StTpcdEdxCorrection::kdZdY); CLRBIT(Mask,StTpcdEdxCorrection::kdXdY); #endif @@ -568,8 +553,6 @@ select firstInnerSectorAnodeWire,lastInnerSectorAnodeWire,numInnerSectorAnodeWir } // sector } // io if (Debug()) Print(); - memset (hist, 0, sizeof(hist)); - memset (checkList, 0, sizeof(checkList)); if (ClusterProfile && GetTFile()) GetTFile()->cd(); #if 0 StMagUtilities::SetDoDistortionT(gFile); @@ -590,29 +573,35 @@ select firstInnerSectorAnodeWire,lastInnerSectorAnodeWire,numInnerSectorAnodeWir {"I","Inner"}, {"O","Outer"} }; - const Name_t PadTime[3] = { - {"Pad","Pad"}, - {"Time","Time"}, - {"Row","Row"}, + const Name_t PadTime[4] = { + {"Pad","Pad - ; Z"}, + {"Time","Time - ; Z"}, + {"PixelPad","Pad ; row ; pad"}, + {"PixelTime","Time ; row ; time "}, }; for (Int_t io = 2; io < 4; io++) { for (Int_t pt = 0; pt < 2; pt++) { TString Name(InOut[io].Name); Name += PadTime[pt].Name; Name += "Mc"; - TString Title(InOut[io].Title); Title += PadTime[pt].Title; Title += "Mc"; + TString Title(InOut[io].Title); Title += PadTime[pt].Title; Title += " Mc"; hist[io][pt] = (TProfile2D *) gDirectory->Get(Name); if (! hist[io][pt]) { - hist[io][pt] = new TProfile2D(Name,Title,nx[pt],xmin[pt],xmax[pt],nz,zmin,zmax,""); - hist[io][pt]->SetMarkerStyle(20); - hist[io][pt]->SetMarkerColor(color++); + hist[io][pt] = new TProfile2D(Name,Title,nx[pt],xmin[pt],xmax[pt],nz,zmin,zmax,""); + hist[io][pt]->SetMarkerStyle(20); + hist[io][pt]->SetMarkerColor(color++); } } } - hist[4][0] = new TProfile2D("dEdxCorSecRow","dEdx correction versus sector and row", - NoOfSectors,0.5,NoOfSectors+0.5, - St_tpcPadConfigC::instance()->numberOfRows(20),0.5,St_tpcPadConfigC::instance()->numberOfRows(20)+0.5,""); - hist[4][1] = new TProfile2D("GainSecRow","Overall gain versus sector and row", - NoOfSectors,0.5,NoOfSectors+0.5, - St_tpcPadConfigC::instance()->numberOfRows(20),0.5,St_tpcPadConfigC::instance()->numberOfRows(20)+0.5,""); + GainHist[0] = new TProfile2D("dEdxCorSecRow","dEdx correction ; sector ; row", + NoOfSectors,0.5,NoOfSectors+0.5, + St_tpcPadConfigC::instance()->numberOfRows(20),0.5,St_tpcPadConfigC::instance()->numberOfRows(20)+0.5,""); + GainHist[1] = new TProfile2D("GainSecRow","Overall gain ; sector ; row", + NoOfSectors,0.5,NoOfSectors+0.5, + St_tpcPadConfigC::instance()->numberOfRows(20),0.5,St_tpcPadConfigC::instance()->numberOfRows(20)+0.5,""); + for (Int_t pt = 0; pt < 2; pt++) { + TString Name(PadTime[pt+2].Name); + TString Title(PadTime[pt+2].Title); Title += " Mc"; + PixelHist[pt] = new TProfile2D(Name,Title, 72, 0.5, 72.5, nx[pt+2],xmin[pt+2],xmax[pt+2],""); + } const Name_t Checks[22] = { {"dEGeant","dE in Geant"}, // 0 {"dSGeant","ds in Geant"}, // 1 @@ -633,9 +622,9 @@ select firstInnerSectorAnodeWire,lastInnerSectorAnodeWire,numInnerSectorAnodeWir {"dS","dS"}, // 16 {"adc","adc"},// 17 {"NE","Total no. of generated electors"}, // 18 - {"dECl","Total log(signal/Nt) in a cluster versus Wire Index"}, // 19 - {"nPdT","log(Total no. of conducting electrons) - log(no. of primary one) versus log(no. primary electrons)"}, // 20 - {"bgVsbg","log10(bg_from_gkin) versus log10(bg_from_mom"} // 21 + {"dECl","Total log(signal/Nt) in a cluster ; Wire Index"}, // 19 + {"nPdT","log(Total no. of conducting electrons) - log(no. of primary one) ; log(no. primary electrons)"}, // 20 + {"bgVsbg","log10(bg_from_gkin) ; log10(bg_from_mom"} // 21 }; const Int_t Npbins = 151; const Int_t NpbinsL = 10; @@ -817,9 +806,6 @@ Int_t StTpcRSMaker::Make(){ // PrintInfo(); continue; } } // special treatment for electron/positron - Int_t qcharge = charge; - if (ipart == 2) qcharge = 101; - if (ipart == 3) qcharge = -101; // Track segment to propagate enum {NoMaxTrackSegmentHits = 100}; static HitPoint_t TrackSegmentHits[NoMaxTrackSegmentHits]; @@ -927,18 +913,22 @@ Int_t StTpcRSMaker::Make(){ // PrintInfo(); checkList[io][3]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),Gain); } // dE/dx correction - Double_t dEdxCor = dEdxCorrection(TrackSegmentHits[iSegHits]); + Double_t dEdxCor = 1.0; + if (TrackSegmentHits[iSegHits].coorLS.position().z() < -0.6) {// prompt hist + dEdxCor = 1.0; + } else { + dEdxCor = dEdxCorrection(TrackSegmentHits[iSegHits]); + } #ifdef __DEBUG__ if (TMath::IsNaN(dEdxCor)) { iBreak++; } #endif - if (dEdxCor <= 0.) continue; + if (dEdxCor < minSignal) continue; if (ClusterProfile) { checkList[io][4]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),dEdxCor); - hist[4][0]->Fill(TrackSegmentHits[iSegHits].Pad.sector(),TrackSegmentHits[iSegHits].Pad.row(),dEdxCor); + GainHist[0]->Fill(TrackSegmentHits[iSegHits].Pad.sector(),TrackSegmentHits[iSegHits].Pad.row(),dEdxCor); } - if (dEdxCor < minSignal) continue; // Initialize propagation Float_t BField[3] = {(Float_t ) TrackSegmentHits[iSegHits].BLS.position().x(), (Float_t ) TrackSegmentHits[iSegHits].BLS.position().y(), @@ -1001,9 +991,9 @@ Int_t StTpcRSMaker::Make(){ // PrintInfo(); Double_t tbkH = TrackSegmentHits[iSegHits].Pad.timeBucket(); tpc_hitC->pad = padH; tpc_hitC->timebucket = tbkH; - pad0 = TMath::Nint(padH + xmin[0]); - tbk0 = TMath::Nint(tbkH + xmin[1]); - Double_t OmegaTau =St_TpcResponseSimulatorC::instance()->OmegaTau()* + pad0 = TMath::Max(1,TMath::Nint(padH + xmin[0])); + tbk0 = TMath::Max(0,TMath::Nint(tbkH + xmin[1])); + Double_t OmegaTau = St_TpcResponseSimulatorC::instance()->OmegaTau()* TrackSegmentHits[iSegHits].BLS.position().z()/5.0;// from diffusion 586 um / 106 um at B = 0/ 5kG Double_t NP = TMath::Abs(tpc_hitC->de)/(St_TpcResponseSimulatorC::instance()->W()*eV* St_TpcResponseSimulatorC::instance()->Cluster()); // from GEANT @@ -1038,14 +1028,7 @@ Int_t StTpcRSMaker::Make(){ // PrintInfo(); mGainLocal = Gain/dEdxCor/NoElPerAdc; // Account dE/dx calibration // end of dE/dx correction // generate electrons: No. of primary clusters per cm - if (mdNdx || mdNdxL10) { - NP = GetNoPrimaryClusters(betaGamma,qcharge); // per cm -#ifdef __DEBUG__ - if (NP <= 0.0) { - iBreak++; continue; - } -#endif - } + NP = StdEdxModel::instance()->dNdx(betaGamma,charge); // per cm if (ClusterProfile) { checkList[io][7]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),NP); } @@ -1057,12 +1040,9 @@ Int_t StTpcRSMaker::Make(){ // PrintInfo(); Float_t dS = 0; Float_t dE = 0; NP = tpc_hitC->dNdx; - static Double_t cLog10 = TMath::Log(10.); if (charge) { dS = zIntDr/NP; - if (mdNdEL10) dE = TMath::Exp(cLog10*mdNdEL10->GetRandom()); - else dE = St_TpcResponseSimulatorC::instance()->W()* - gRandom->Poisson(St_TpcResponseSimulatorC::instance()->Cluster()); + dE = StdEdxModel::instance()->dNdE(); } else { // charge == 0 geantino // for Laserino assume dE/dx = 25 keV/cm; dE = 10; // eV @@ -1121,8 +1101,6 @@ Int_t StTpcRSMaker::Make(){ // PrintInfo(); TotalSignalInCluster = 0; Int_t WireIndex = 0; for (Int_t ie = 0; ie < Nt; ie++) { - tpc_hitC->ne++; - QAv = mPolya[io]->GetRandom(); // transport to wire gRandom->Rannor(rX,rY); StTpcLocalSectorCoordinate xyzE(xyzC.x()+rX*sigmaT, @@ -1141,6 +1119,19 @@ Int_t StTpcRSMaker::Make(){ // PrintInfo(); #endif TCL::vadd(xyzE.position().xyz(),xyzR.GetArray(),xyzE.position().xyz(),3); } + // Check Gatting Grid + Double_t zGG = xyzE.position().z(); + Double_t GGTransperency = 1.0; + if (zGG > -0.6) { // non Prompt hits before Cathode wire plane + Double_t driftTime = 1e6*zGG/driftVelocity; // usec + driftTime -= St_starTriggerDelayC::instance()->TrigT0GG(io); // trigger delay + Gating Grid cables + Double_t lGGTransperency = St_GatingGridBC::instance()->CalcCorrection(0,driftTime); + if (lGGTransperency < -9) continue; + GGTransperency = TMath::Exp(lGGTransperency); + if (GGTransperency < minSignal) continue; + } + tpc_hitC->ne++; + QAv = GGTransperency*mPolya[io]->GetRandom(); Double_t y = xyzE.position().y(); Double_t alphaVariation = InnerAlphaVariation[sector-1]; // Transport to wire @@ -1219,12 +1210,16 @@ Int_t StTpcRSMaker::Make(){ // PrintInfo(); if (ClusterProfile) { if (TotalSignal > 0) { if (hist[ioH][0]) { - for (Int_t p = 0; p < kPadMax; p++) + for (Int_t p = 0; p < kPadMax; p++) { hist[ioH][0]->Fill((p+pad0)-padH,TrackSegmentHits[iSegHits].xyzG.position().z(),padsdE[p]/TotalSignal); + PixelHist[0]->Fill(row, p+pad0,padsdE[p]/TotalSignal); + } } if (hist[ioH][1]) { - for (Int_t t = 0; t < kTimeBacketMax; t++) + for (Int_t t = 0; t < kTimeBacketMax; t++) { hist[ioH][1]->Fill((t+tbk0+0.5)-tbkH,TrackSegmentHits[iSegHits].xyzG.position().z(),tbksdE[t]/TotalSignal); + PixelHist[1]->Fill(row, (t+tbk0+0.5),tbksdE[t]/TotalSignal); + } } } checkList[io][15]->Fill(TrackSegmentHits[iSegHits].xyzG.position().z(),tpc_hitC->de); @@ -1277,37 +1272,6 @@ Int_t StTpcRSMaker::Make(){ // PrintInfo(); return kStOK; } //________________________________________________________________________________ -Double_t StTpcRSMaker::GetNoPrimaryClusters(Double_t betaGamma, Int_t charge) { - Double_t beta = betaGamma/TMath::Sqrt(1.0 + betaGamma*betaGamma); -#if defined(Old_dNdx_Table) - Double_t dNdx = 1.21773e+01*Bichsel::Instance()->GetI70M(TMath::Log10(betaGamma)); -#else -#if defined(ElectronHack) - Int_t elepos = charge/100; - Double_t dNdx = 0; - if (mdNdx ) dNdx = mdNdx->Interpolate(betaGamma); - else if (mdNdxL10) dNdx = mdNdxL10->Interpolate(TMath::Log10(betaGamma)); - if (elepos) { - dNdx += 1.21773e+01*Bichsel::Instance()->GetI70M(TMath::Log10(betaGamma)); - dNdx /= 2; - } -#else /* new H.Bichsel dNdx table 09/12/11 */ - Double_t dNdx = 0; - if (mdNdx ) dNdx = mdNdx->Interpolate(betaGamma); - else if (mdNdxL10) dNdx = mdNdxL10->Interpolate(TMath::Log10(betaGamma)); -#endif /* Old_dNdx_Table || ElectronHack */ -#endif - Double_t Q_eff = TMath::Abs(charge%100); - if (Q_eff > 1) { - // Effective charge from GEANT gthion.F - Double_t w1 = 1.034 - 0.1777*TMath::Exp(-0.08114*Q_eff); - Double_t w2 = beta*TMath::Power(Q_eff,-2./3.); - Double_t w3 = 121.4139*w2 + 0.0378*TMath::Sin(190.7165*w2); - Q_eff *= 1. -w1*TMath::Exp(-w3); - } - return Q_eff*Q_eff*dNdx; -} -//________________________________________________________________________________ Double_t StTpcRSMaker::ShaperFunc(Double_t *x, Double_t *par) { Double_t tau = par[0]; Double_t width = par[1]; @@ -2059,7 +2023,7 @@ void StTpcRSMaker::GenerateSignal(HitPoint_t &TrackSegmentHits, Int_t sector, In } if (ClusterProfile) { checkList[io][12]->Fill(TrackSegmentHits.xyzG.position().z(),gain); - hist[4][1]->Fill(sector,row,gain); + GainHist[1]->Fill(sector,row,gain); } // Double_t localXDirectionCoupling = localXDirectionCouplings[pad-padMin]; Double_t localXDirectionCoupling = gain*XDirectionCouplings[pad-padMin]; @@ -2181,324 +2145,3 @@ Double_t StTpcRSMaker::dEdxCorrection(HitPoint_t &TrackSegmentHits) { } //________________________________________________________________________________ #undef PrPP -//________________________________________________________________________________ -// $Id: StTpcRSMaker.cxx,v 1.92 2020/05/22 20:49:19 fisyak Exp $ -// $Log: StTpcRSMaker.cxx,v $ -// Revision 1.92 2020/05/22 20:49:19 fisyak -// Wrong alarm, take it back -// -// Revision 1.91 2020/05/22 20:15:15 fisyak -// Fix bug in fgTriggerT0 (seconds <=> microseconds) -// -// Revision 1.90 2020/05/17 15:43:49 fisyak -// Acount StTpcBXT0CorrEPDC correction -// -// Revision 1.89 2019/05/22 21:30:58 fisyak -// Fix bug3390 (thanks to Irakli), add St_TpcAdcCorrectionMDF -// -// Revision 1.88 2019/04/29 20:11:21 fisyak -// Fix for TrackDirection, add extra correction for the 1st pad row -// -// Revision 1.87 2019/04/18 14:02:23 fisyak -// Keep hits with bad dE/dx information, revisit handling of pedRMS for Tpx and iTPC -// -// Revision 1.86 2018/12/16 14:26:30 fisyak -// Switch off DEBUG, use local position for phi correction -// -// Revision 1.85 2018/12/09 23:22:59 fisyak -// Reshape -// -// Revision 1.84 2018/11/29 22:19:49 fisyak -// Restore __STOPPED_ELECTRONS__, split for Inner and Outer sectors, adjusted gain for Run XVIII -// -// Revision 1.83 2018/11/20 19:51:15 fisyak -// Temporarely disable __STOPPED_ELECTRONS__ to check effect of this on no. of primary tracks for 2010 AuAu200 sample -// -// Revision 1.82 2018/11/05 01:05:19 fisyak -// Replace assert to error message -// -// Revision 1.81 2018/11/01 13:27:20 fisyak -// Synchronize mCutEle with GEANT3 tracking media setting for TPCE_SENSITIVE_GAS, bug#3369 -// -// Revision 1.80 2018/10/17 20:45:28 fisyak -// Restore update for Run XVIII dE/dx calibration removed by Gene on 08/07/2018 -// -// Revision 1.77 2018/02/18 21:08:54 perev -// Move back DSmirnov correction -// -// Revision 1.75 2017/02/14 23:40:35 fisyak -// Add new Table to correct dE/dx pad dependence, 2017 dAu20-62 calibration -// -// Revision 1.74 2016/12/29 16:30:56 fisyak -// make switch to account __STOPPED_ELECTRONS__ -// -// Revision 1.73 2016/12/19 15:22:39 fisyak -// Fix typo -// -// Revision 1.72 2016/12/19 15:20:20 fisyak -// Fix bug 3268: add missing correction for TpcAvgPowerSupply, add treatment for stopping electrons and gammas -// -// Revision 1.71 2016/09/18 22:45:25 fisyak -// Clean up, add Heed model, adjust for new StTpcdEdxCorrections -// -// Revision 1.71 2015/10/08 13:58:45 fisyak -// Add requirement for Check Plots for TTree file -// -// Revision 1.70 2015/07/19 22:14:07 fisyak -// Clean up __PAD_BLOCK__, recalculate no. of real hits in g2t_track n_tpc_hitC (excluding pseudo pad row), add current and accumulated charge in dE/dx correction -// -// Revision 1.69 2014/10/21 15:33:48 fisyak -// Clean up, fix bug found by gcc482 -// -// Revision 1.68 2014/07/27 13:26:26 fisyak -// Add cast for c++11 option -// -// Revision 1.67 2013/02/01 15:54:09 fisyak -// Add handle for separate Inner and Outer sector time off set -// -// Revision 1.66 2012/11/13 20:46:16 fisyak -// Add wider Voltage range for accepted clusteds (-500V) than for dEdx calculation -// -// Revision 1.65 2012/10/23 20:08:57 fisyak -// Add corrections for iTpx upgrade -// -// Revision 1.64 2012/09/27 19:17:02 fisyak -// Fix missing declaration -// -// Revision 1.63 2012/09/27 16:14:43 fisyak -// Change debug print out scheme -// -// Revision 1.62 2012/09/13 21:54:43 fisyak -// replace elsif by else and if -// -// Revision 1.61 2012/09/13 21:02:52 fisyak -// Corrections for iTpx -// -// Revision 1.60 2012/06/04 15:14:18 fisyak -// restore old hack for dN/dx table to fix bug #2347 -// -// Revision 1.59 2012/05/07 15:36:22 fisyak -// Remove hardcoded TPC parameters -// -// Revision 1.58 2012/04/03 14:05:18 fisyak -// Speed up using GetSaveL (__PAD_BLOCK__), sluggish shape histograms, Heed electron generation -// -// Revision 1.57 2011/12/23 16:38:25 fisyak -// Remove __DEBUG__ and __ClusterProfile__ from default, reduce arrays and add check for bounds -// -// Revision 1.55 2011/12/20 21:09:56 fisyak -// change defaults: shark measurements: old default => 46.6%, wire histograms => 38.9%, wire map => 12.5 + 10.2, pad block => 15% -// -// Revision 1.54 2011/12/13 17:23:22 fisyak -// remove YXTProd, add WIREHISTOGRAM and WIREMAP, use particle definition from StarClassLibrary -// -// Revision 1.53 2011/10/14 23:27:51 fisyak -// Back to standard version -// -// Revision 1.51 2011/09/21 15:34:31 fisyak -// Restore K3IP parameter -// -// Revision 1.50 2011/09/18 22:39:48 fisyak -// Extend dN/dx table (H.Bichsel 09/12/2011) to fix bug #2174 and #2181, clean-up -// -// Revision 1.49 2011/04/05 20:55:02 fisyak -// Fix betaMin calculations (thanx VP) -// -// Revision 1.48 2011/03/17 14:29:31 fisyak -// Add extrapolation in region beta*gamma < 0.3 -// -// Revision 1.47 2011/03/13 15:46:44 fisyak -// Replace step function by interpolation for dN/dx versus beta*gamma -// -// Revision 1.46 2011/02/23 20:14:31 perev -// Hack to avoid sqrt(-) -// -// Revision 1.45 2010/12/16 15:36:07 fisyak -// cut hits outside time buckets range -// -// Revision 1.44 2010/12/01 20:59:46 fisyak -// Remove special treatment for delta-electrons, this will cause that IdTruth for cluster will be degradated because charge from delta-electrons will be accounted with delta-electons track Id but not with original particle Id as was before -// -// Revision 1.43 2010/10/28 23:42:34 fisyak -// extra t0 off set for Altro chip -// -// Revision 1.42 2010/10/22 18:13:33 fisyak -// Add fix from Lokesh AuAu7 2010 embdedding -// -// Revision 1.41 2010/09/01 23:12:01 fisyak -// take out __ClusterProfile__ -// -// Revision 1.40 2010/06/14 23:34:26 fisyak -// Freeze at Version V -// -// Revision 1.39 2010/05/24 16:11:03 fisyak -// Return back to time simulation for each pad, organize parameters into TpcResponseSimulator table -// -// Revision 1.38 2010/04/24 19:58:54 fisyak -// swap shift sign -// -// Revision 1.37 2010/04/24 15:56:32 fisyak -// Jan found shift in z by one time bucket -// -// Revision 1.36 2010/04/20 13:56:24 fisyak -// Switch off __ClusterProfile__ -// -// Revision 1.35 2010/04/16 19:29:35 fisyak -// W is in eV now -// -// Revision 1.34 2010/04/01 22:17:06 fisyak -// Add checking for TPC is switched off at all and stop if so -// -// Revision 1.33 2010/03/22 23:45:05 fisyak -// Freeze version with new parameters table -// -// Revision 1.32 2010/03/16 19:41:46 fisyak -// Move diffusion and sec/row correction in DB, clean up -// -// Revision 1.31 2010/03/02 21:10:27 fisyak -// Make aware about TpcRDOMasks -// -// Revision 1.30 2010/02/26 18:53:33 fisyak -// Take longitudinal Diffusion from Laser track fit, add Gating Grid -// -// Revision 1.29 2010/02/16 00:21:23 fisyak -// Speed up by a factor 3.5 by ignoring individual pad T0 -// -// Revision 1.28 2010/01/26 19:47:25 fisyak -// Include dE/dx calibration and distortions in the simulation -// -// Revision 1.27 2009/11/25 21:32:52 fisyak -// Comment out cluster profile histograms -// -// Revision 1.26 2009/11/03 22:38:53 fisyak -// Freeze version rcf9108.J -// -// Revision 1.25 2009/10/30 21:12:00 fisyak -// Freeze version rcf9108.F, Clean up -// -// Revision 1.24 2009/10/26 18:50:58 fisyak -// Clean up from Bichel's stuff -// -// Revision 1.23 2009/10/12 23:54:12 fisyak -// Restore T0Jitter, remove differential in Tpx signal -// -// Revision 1.22 2009/10/03 21:29:09 fisyak -// Clean up, move all TpcT related macro into StTpcMcAnalysisMaker -// -// Revision 1.21 2009/10/01 14:53:06 fisyak -// Add T0Jitter -// -// Revision 1.20 2009/09/27 01:30:48 fisyak -// Restate T0Jitter -// -// Revision 1.19 2009/09/27 01:24:58 fisyak -// Restate T0Jitter -// -// Revision 1.18 2009/09/21 13:20:39 fisyak -// Variant O4, no mSigmaJitter, 100 keV -// -// Revision 1.17 2009/09/01 15:06:44 fisyak -// Version N -// -// Revision 1.16 2009/08/25 20:39:40 fisyak -// Variant K -// -// Revision 1.15 2009/08/25 15:45:58 fisyak -// Version J -// -// Revision 1.14 2009/08/24 20:16:41 fisyak -// Freeze with new Altro parameters -// -// Revision 1.13 2008/12/29 15:24:54 fisyak -// Freeze ~/WWW/star/Tpc/TpcRS/ComparisonMIP31 -// -// Revision 1.12 2008/12/18 23:06:37 fisyak -// Take care about references to TGiant -// -// Revision 1.11 2008/12/12 21:41:41 fisyak -// Freeze -// -// Revision 1.10 2008/10/06 19:10:23 fisyak -// BichlePPMIP3 -// -// Revision 1.9 2008/10/03 20:25:29 fisyak -// Version BichselMIP2 -// -// Revision 1.8 2008/08/19 16:01:15 fisyak -// Version 21 -// -// Revision 1.7 2008/08/18 15:54:25 fisyak -// Version 20 -// -// Revision 1.6 2008/07/30 23:53:19 fisyak -// Freeze -// -// Revision 1.5 2008/07/18 16:22:50 fisyak -// put a factor 2.5 for tauIntegration -// -// Revision 1.3 2008/06/25 20:02:32 fisyak -// The first set of parametrs for Altro, Remove gains for the moment -// -// Revision 1.2 2008/06/19 22:45:43 fisyak -// Freeze problem with TPX parameterization -// -// Revision 1.1.1.1 2008/04/28 14:39:47 fisyak -// Start new Tpc Response Simulator -// -// Revision 1.20 2008/04/24 10:42:03 fisyak -// Fix binning issues -// -// Revision 1.19 2008/04/04 15:00:11 fisyak -// Freeze before shaper modifications -// -// Revision 1.18 2005/02/07 21:40:31 fisyak -// rename antique TGeant3 to TGiant3 -// -// Revision 1.17 2005/01/26 23:28:38 fisyak -// Check boundary for sorted tpc_hit array -// -// Revision 1.16 2005/01/26 21:45:31 fisyak -// Freeze correction made in June -// -// Revision 1.14 2004/06/04 17:09:02 fisyak -// Change tau in Chaper and OmegaTau for gas -// -// Revision 1.13 2004/05/29 21:16:27 fisyak -// Fix pad direction, add sorting for ADC/cluster nonlinearity, replace product by sum of logs -// -// Revision 1.12 2004/05/17 19:45:08 fisyak -// Clean up, add pseudo padrows -// -// Revision 1.11 2004/05/05 17:41:52 fisyak -// Take K3 from E.Mathieson book -// -// Revision 1.10 2004/05/04 13:39:06 fisyak -// Add TF1 -// -// Revision 1.9 2004/05/02 20:54:18 fisyak -// fix t0 offset -// -// Revision 1.8 2004/04/22 01:05:03 fisyak -// Freeze the version before modification parametrization for K3 -// -// Revision 1.7 2004/04/12 14:30:07 fisyak -// Propagate cluster as a whole -// -// Revision 1.6 2004/04/06 01:50:13 fisyak -// Switch from Double_t to Float_t for sum -// -// Revision 1.5 2004/03/30 19:30:04 fisyak -// Add Laser -// -// Revision 1.4 2004/03/21 19:00:43 fisyak -// switch to GEANT step length -// -// Revision 1.3 2004/03/20 17:57:15 fisyak -// Freeze the version of PAI model -// -// Revision 1.2 2004/03/17 20:47:43 fisyak -// Add version with TrsCluster TTree -// -// Revision 1.1.1.1 2004/03/05 20:51:25 fisyak -// replacement for Trs -// diff --git a/StRoot/StTpcRSMaker/StTpcRSMaker.h b/StRoot/StTpcRSMaker/StTpcRSMaker.h index 66e3bbf6c60..ada9282164b 100644 --- a/StRoot/StTpcRSMaker/StTpcRSMaker.h +++ b/StRoot/StTpcRSMaker/StTpcRSMaker.h @@ -58,7 +58,6 @@ class StTpcRSMaker : public StMaker { TF1 *GetTimeShape0(Int_t io = 0) {return fgTimeShape0[io];} TF1 *GetTimeShape3(Int_t io = 0) {return fgTimeShape3[io];} TF1 *GetHeed() {return mHeed;} - Double_t GetNoPrimaryClusters(Double_t betaGamma, Int_t charge); virtual void Print(Option_t *option="") const; StTpcDigitalSector *DigitizeSector(Int_t sector); void SetLaserScale(Double_t m=1) {mLaserScale = m;} @@ -100,9 +99,6 @@ class StTpcRSMaker : public StMaker { Char_t beg[1]; //! TTree *fTree; //! SignalSum_t *m_SignalSum; //! - TH1D* mdNdx; //! - TH1D* mdNdxL10; //! - TH1D* mdNdEL10; //! TF1F *mShaperResponses[2][24]; //! TF1F *mChargeFraction[2][24]; //! TF1F *mPadResponseFunction[2][24]; //! diff --git a/StRoot/StdEdxY2Maker/StTpcdEdxCorrection.cxx b/StRoot/StdEdxY2Maker/StTpcdEdxCorrection.cxx index 54a688beab5..f65e1007d8d 100644 --- a/StRoot/StdEdxY2Maker/StTpcdEdxCorrection.cxx +++ b/StRoot/StdEdxY2Maker/StTpcdEdxCorrection.cxx @@ -31,8 +31,6 @@ #include "StDetectorDbMaker/St_tpcPressureBC.h" #include "StDetectorDbMaker/St_TpcDriftDistOxygenC.h" #include "StDetectorDbMaker/St_TpcMultiplicityC.h" -#include "StDetectorDbMaker/St_GatingGridC.h" -#include "StDetectorDbMaker/St_GatingGridBC.h" #include "StDetectorDbMaker/St_TpcZCorrectionBC.h" #include "StDetectorDbMaker/St_TpcZCorrectionCC.h" #include "StDetectorDbMaker/St_tpcMethaneInC.h" @@ -138,7 +136,6 @@ void StTpcdEdxCorrection::ReSetCorrections() { assert(tpcGas); } SettpcGas(tpcGas); - memset (m_Corrections, 0, sizeof(m_Corrections)); mTimeBinWidth = 1./StTpcDb::instance()->Electronics()->samplingFrequency(); mInnerSectorzOffset = StTpcDb::instance()->Dimensions()->zInnerOffset(); mOuterSectorzOffset = StTpcDb::instance()->Dimensions()->zOuterOffset(); @@ -157,28 +154,37 @@ void StTpcdEdxCorrection::ReSetCorrections() { m_Corrections[kAdcCorrectionC ] = dEdxCorrection_t("TpcAdcCorrectionC" ,"alternative ADC/Clustering nonlinearity correction" ,St_TpcAdcCorrectionCC::instance()); m_Corrections[kEdge ] = dEdxCorrection_t("TpcEdge" ,"Gain on distance from Chamber edge" ,St_TpcEdgeC::instance()); m_Corrections[kAdcCorrectionMDF ] = dEdxCorrection_t("TpcAdcCorrectionMDF" ,"ADC/Clustering nonlinearity correction MDF" ,St_TpcAdcCorrectionMDF::instance()); + /* used only for tests m_Corrections[kAdcCorrection3MDF ] = dEdxCorrection_t("TpcAdcCorrection3MDF","ADC/Clustering nonlinearity correction MDF 3D" ,St_TpcAdcCorrection3MDF::instance()); m_Corrections[kAdcCorrection4MDF ] = dEdxCorrection_t("TpcAdcCorrection4MDF","ADC/Clustering nonlinearity correction MDF 4D" ,St_TpcAdcCorrection4MDF::instance()); m_Corrections[kAdcCorrection5MDF ] = dEdxCorrection_t("TpcAdcCorrection5MDF","ADC/Clustering nonlinearity correction MDF+4D" ,St_TpcAdcCorrection5MDF::instance()); + */ m_Corrections[kAdcCorrection6MDF ] = dEdxCorrection_t("TpcAdcCorrection6MDF","alternative ADC/Clustering nonlinearity correction MDF+4D" ,St_TpcAdcCorrection6MDF::instance()); m_Corrections[kTpcdCharge ] = dEdxCorrection_t("TpcdCharge" ,"ADC/Clustering undershoot correction" ,St_TpcdChargeC::instance()); + /* m_Corrections[kTpcrCharge ] = dEdxCorrection_t("TpcrCharge" ,"ADC/Clustering rounding correction" ,St_TpcrChargeC::instance()); + */ m_Corrections[kTpcCurrentCorrection ] = dEdxCorrection_t("TpcCurrentCorrection","Correction due to sagg of Voltage due to anode current" ,St_TpcCurrentCorrectionC::instance()); m_Corrections[kTpcRowQ ] = dEdxCorrection_t("TpcRowQ" ,"Gas gain correction for row versus accumulated charge," ,St_TpcRowQC::instance()); m_Corrections[kTpcAccumulatedQ ] = dEdxCorrection_t("TpcAccumulatedQ" ,"Gas gain correction for HV channel versus accumulated charge," ,St_TpcAccumulatedQC::instance()); m_Corrections[kTpcSecRowB ] = dEdxCorrection_t("TpcSecRowB" ,"Gas gain correction for sector/row" ,St_TpcSecRowBC::instance()); + /* m_Corrections[kTpcSecRowC ] = dEdxCorrection_t("TpcSecRowC" ,"Additional Gas gain correction for sector/row" ,St_TpcSecRowCC::instance()); + */ m_Corrections[ktpcPressure ] = dEdxCorrection_t("tpcPressureB" ,"Gain on Gas Density due to Pressure" ,St_tpcPressureBC::instance()); m_Corrections[ktpcTime ] = dEdxCorrection_t("tpcTime" ,"Unregognized time dependce" ,St_tpcTimeDependenceC::instance()); m_Corrections[kDrift ] = dEdxCorrection_t("TpcDriftDistOxygen" ,"Correction for Electron Attachment due to O2" ,St_TpcDriftDistOxygenC::instance()); m_Corrections[kMultiplicity ] = dEdxCorrection_t("TpcMultiplicity" ,"Global track multiplicity dependence" ,St_TpcMultiplicityC::instance()); - m_Corrections[kGatingGrid ] = dEdxCorrection_t("GatingGrid" ,"Variation due to Gating Grid transperancy" ,St_GatingGridBC::instance()); m_Corrections[kzCorrectionC ] = dEdxCorrection_t("TpcZCorrectionC" ,"Variation on drift distance with Gating Grid one" ,St_TpcZCorrectionCC::instance()); m_Corrections[kzCorrection ] = dEdxCorrection_t("TpcZCorrectionB" ,"Variation on drift distance without Gating Gird one" ,St_TpcZCorrectionBC::instance()); + /* Deactivated m_Corrections[ktpcMethaneIn ] = dEdxCorrection_t("tpcMethaneIn" ,"Gain on Methane content" ,St_tpcMethaneInC::instance()); - m_Corrections[ktpcGasTemperature ] = dEdxCorrection_t("tpcGasTemperature" ,"Gain on Gas Dens. due to Temperature" ,St_tpcGasTemperatureC::instance()); + */ + m_Corrections[ktpcGasTemperature ] = dEdxCorrection_t("tpcGasTemperature" ,"Gain on Gas Dens. due to Output (T5) Temperature" ,St_tpcGasTemperatureC::instance()); + /* Deactivated m_Corrections[ktpcWaterOut ] = dEdxCorrection_t("tpcWaterOut" ,"Gain on Water content" ,St_tpcWaterOutC::instance()); m_Corrections[kSpaceCharge ] = dEdxCorrection_t("TpcSpaceCharge" ,"Gain on space charge near the wire" ,St_TpcSpaceChargeC::instance()); + */ m_Corrections[kPhiDirection ] = dEdxCorrection_t("TpcPhiDirection" ,"Gain on interception angle" ,St_TpcPhiDirectionC::instance()); m_Corrections[kTanL ] = dEdxCorrection_t("TpcTanL" ,"Gain on Tan(lambda)" ,St_TpcTanLC::instance()); m_Corrections[kdXCorrection ] = dEdxCorrection_t("TpcdXCorrectionB" ,"dX correction" ,St_TpcdXCorrectionBC::instance()); @@ -189,11 +195,15 @@ void StTpcdEdxCorrection::ReSetCorrections() { m_Corrections[kTpcZDC ] = dEdxCorrection_t("TpcZDC" ,"Gain on Zdc CoincidenceRate" ,St_TpcZDCC::instance()); m_Corrections[kTpcPadMDF ] = dEdxCorrection_t("TpcPadCorrectionMDF" ,"Gain Variation along the anode wire" ,St_TpcPadCorrectionMDF::instance()); m_Corrections[kTpcPadMDC ] = dEdxCorrection_t("TpcPadCorrectionMDC" ,"Gain Variation along the anode wire with track curvature" ,St_TpcPadCorrectionMDC::instance()); + /* m_Corrections[kAdcI ] = dEdxCorrection_t("TpcAdcI" ,"Gain on Accumulated Adc on a socket)" ,St_TpcAdcIC::instance()); + */ m_Corrections[knPad ] = dEdxCorrection_t("TpcnPad" ,"Gain on cluster length in pads" ,St_TpcnPadC::instance()); - m_Corrections[knTbk ] = dEdxCorrection_t("TpcnTbk" ,"Gain on cluster length i time buckets" ,St_TpcnTbkC::instance()); + m_Corrections[knTbk ] = dEdxCorrection_t("TpcnTbk" ,"Gain on cluster length in time buckets" ,St_TpcnTbkC::instance()); + /* m_Corrections[kdZdY ] = dEdxCorrection_t("TpcdZdY" ,"Gain on track dZ/dY" ,St_TpcdZdYC::instance()); m_Corrections[kdXdY ] = dEdxCorrection_t("TpcdXdY" ,"Gain on track dX/dY" ,St_TpcdXdYC::instance()); + */ m_Corrections[kTpcLast ] = dEdxCorrection_t("Final" ,"" ,0); m_Corrections[kTpcLengthCorrection ] = dEdxCorrection_t("TpcLengthCorrectionB" ,"Variation vs Track length and relative error in Ionization" ,St_TpcLengthCorrectionBC::instance()); m_Corrections[kTpcLengthCorrectionMDF] = dEdxCorrection_t("TpcLengthCorrectionMDF","Variation vs Track length and and rel. error in dE/dx" ,St_TpcLengthCorrectionMDF::instance()); @@ -262,7 +272,7 @@ void StTpcdEdxCorrection::ReSetCorrections() { } // Check consistency of active chairs if ( m_Corrections[kzCorrectionC].Chair) { // if kzCorrectionC is already active - vector kvect = {kzCorrection, kGatingGrid}; + vector kvect = {kzCorrection}; // , kGatingGrid}; for (auto k : kvect) { if (m_Corrections[k].Chair) { LOG_WARN << m_Corrections[kzCorrectionC].Name << " is already active => disable " << m_Corrections[k].Name << endm; @@ -388,19 +398,16 @@ Int_t StTpcdEdxCorrection::dEdxCorrection(dEdxY2_t &CdEdx, Bool_t doIT) { #endif #endif Double_t ZdriftDistance = CdEdx.ZdriftDistance; - Double_t driftDistance2GG = ZdriftDistance; ESector kTpcOutIn = kTpcOuter; Float_t gasGain = 1; Float_t gainNominal = 0; St_tss_tssparC *tsspar = St_tss_tssparC::instance(); if (row <= St_tpcPadConfigC::instance()->innerPadRows(sector)) { kTpcOutIn = kTpcInner; - driftDistance2GG += mInnerSectorzOffset; gainNominal = tsspar->gain_in()*tsspar->wire_coupling_in(); gasGain = tsspar->gain_in(sector,row) *tsspar->wire_coupling_in(); } else { kTpcOutIn = kTpcOuter; - driftDistance2GG += mOuterSectorzOffset; gainNominal = tsspar->gain_out()*tsspar->wire_coupling_out(); gasGain = tsspar->gain_out(sector,row)*tsspar->wire_coupling_out(); } @@ -409,7 +416,7 @@ Int_t StTpcdEdxCorrection::dEdxCorrection(dEdxY2_t &CdEdx, Bool_t doIT) { iok = 1; return iok; } - CdEdx.driftTime = driftDistance2GG/gStTpcDb->DriftVelocity(sector)*1e6 - m_TrigT0;// musec + CdEdx.driftTime = ZdriftDistance/gStTpcDb->DriftVelocity(sector)*1e6 - m_TrigT0;// musec // Double_t gainAVcorr = gasGain/gainNominal; mAdc2GeV = tsspar->ave_ion_pot() * tsspar->scale()/gainNominal; Double_t Adc2GeVReal = tsspar->ave_ion_pot() * tsspar->scale()/gasGain; @@ -469,7 +476,7 @@ Int_t StTpcdEdxCorrection::dEdxCorrection(dEdxY2_t &CdEdx, Bool_t doIT) { if (! table) goto ENDL; const TpcSecRowCor_st *gain = table->GetTable() + sector - 1; gc = gain->GainScale[row-1]; - // gcRMS = gain->GainRms[row-1]; + // gcRMS = gain->GainRms[row-1];1 // if (gc <= 0.0 || gcRMS <= 0.0) { if (gc <= 0.0) { return k; @@ -542,18 +549,6 @@ Int_t StTpcdEdxCorrection::dEdxCorrection(dEdxY2_t &CdEdx, Bool_t doIT) { } else if (k == kAdcCorrection6MDF) { goto ENDL; // kAdcCorrection6MDF is in kAdcCorrectionC } -#if 0 - if (k == kGatingGrid) { - // Take care about prompt hits and Gating Grid region in Simulation - if (ZdriftDistance <= 0.0) goto ENDL; // prompt hits - Double_t dEcor = ((St_GatingGridC *)m_Corrections[k].Chair)->CalcCorrection(l,VarXs[kGatingGrid]); - if (dEcor < -9.9) { - return k; - } - dE *= TMath::Exp(-dEcor); - goto ENDL; - } -#endif cor = ((St_tpcCorrection *) m_Corrections[k].Chair->Table())->GetTable(); if (! cor) goto ENDL; NLoops = cor->type/100000 + 1; @@ -636,8 +631,11 @@ Int_t StTpcdEdxCorrection::dEdxCorrection(dEdxY2_t &CdEdx, Bool_t doIT) { return k; } } - if (VarXs[k] < 0) {// prompt hits + if (ZdriftDistance < -0.6) {// prompt hits after cathode wire plane dE *= TMath::Exp(1.2); + goto ENDL; + } else if (VarXs[k] < 0.0) {// hits after Gating Grid and before cathode wire plane + return k; } if (k == kzCorrectionC && corl->type == 20) { Int_t np = TMath::Abs(corl->npar)%100; @@ -671,27 +669,10 @@ Int_t StTpcdEdxCorrection::dEdxCorrection(dEdxY2_t &CdEdx, Bool_t doIT) { return k; } } - } - if (corl->npar%100) { - Double_t dECor = 0; - if (k == kGatingGrid) { - dECor = TMath::Exp(-((St_GatingGridBC *)m_Corrections[k].Chair)->CalcCorrection(l,VarXs[k])); - } else { - dECor = TMath::Exp(-chairC->CalcCorrection(l,VarXs[k])); - } -#if 0 - if (TMath::IsNaN(dECor)) { - static Int_t iBreak = 0; - iBreak++; + if (corl->npar%100) { + Double_t dECor = TMath::Exp(-chairC->CalcCorrection(l,VarXs[k])); + dE *= dECor; } -#endif - dE *= dECor; -#if 0 - if (TMath::IsNaN(dE)) { - static Int_t iBreak = 0; - iBreak++; - } -#endif } else { // repeatable for (m = 0; m < NLoops; m++, l += nrows) { corl = cor + l; @@ -752,6 +733,7 @@ Int_t StTpcdEdxCorrection::dEdxTrackCorrection(EOptions opt, Int_t type, dst_ded if (! m_Corrections[k].Chair) return 0; Int_t l = 2*type; Int_t nrows = 0; + Double_t dEdxError = 0; const St_tpcCorrectionC *chairC = dynamic_cast(m_Corrections[k].Chair); switch (k) { case kTpcLengthCorrection: @@ -762,7 +744,8 @@ Int_t StTpcdEdxCorrection::dEdxTrackCorrection(EOptions opt, Int_t type, dst_ded case 2: // fit if (nrows > l+1) { dedx.dedx[0] *= TMath::Exp(-chairC->CalcCorrection( l,LogTrackLength)); - dedx.dedx[1] = chairC->CalcCorrection(l+1,LogTrackLength); + dEdxError = chairC->CalcCorrection(l+1,LogTrackLength); + if (dEdxError > 0) dedx.dedx[1] = dEdxError; } if (nrows > l+6) { dedx.dedx[0] *= TMath::Exp(-chairC->CalcCorrection(l+6,LogTrackLength)); @@ -783,7 +766,8 @@ Int_t StTpcdEdxCorrection::dEdxTrackCorrection(EOptions opt, Int_t type, dst_ded case 2: // fit if (nrows > l+1) { dedx.dedx[0] *= TMath::Exp(-((St_MDFCorrectionC *)m_Corrections[k].Chair)->Eval( l,xx)); - dedx.dedx[1] = ((St_MDFCorrectionC *)m_Corrections[k].Chair)->Eval(l+1,xx); + dEdxError = ((St_MDFCorrectionC *)m_Corrections[k].Chair)->Eval(l+1,xx); + if (dEdxError > 0) dedx.dedx[1] = dEdxError; } break; default: diff --git a/StRoot/StdEdxY2Maker/StTpcdEdxCorrection.h b/StRoot/StdEdxY2Maker/StTpcdEdxCorrection.h index 66048c4cf59..c2009ba0fa3 100644 --- a/StRoot/StdEdxY2Maker/StTpcdEdxCorrection.h +++ b/StRoot/StdEdxY2Maker/StTpcdEdxCorrection.h @@ -25,7 +25,7 @@ struct dE_t { }; //________________________________________________________________________________ struct dEdxCorrection_t { - dEdxCorrection_t(const Char_t *name = "",const Char_t *title = "", TChair *chair=0, Int_t n=0) : + dEdxCorrection_t(const Char_t *name = "",const Char_t *title = "", TChair *chair=0) : Name(name), Title(title), Chair(chair) {} const Char_t *Name; const Char_t *Title; @@ -62,7 +62,7 @@ class StTpcdEdxCorrection : public TObject { kzCorrection = 22,//Z ktpcMethaneIn = 23,//m ktpcGasTemperature = 24,//T - ktpcWaterOut = 25,//W 7 + ktpcWaterOut = 25,//W kSpaceCharge = 26,//C space charge near the wire kPhiDirection = 27,//p correction wrt local interception angle kTanL = 28,//p correction wrt local tan(lambda) diff --git a/StarDb/Calibrations/tpc/GatingGridB.20190225.230040.C b/StarDb/Calibrations/tpc/GatingGridB.20190225.230040.C deleted file mode 100644 index 4240f0a88d3..00000000000 --- a/StarDb/Calibrations/tpc/GatingGridB.20190225.230040.C +++ /dev/null @@ -1,57 +0,0 @@ -#include "tables/St_tpcCorrection_Table.h" - -TDataSet *CreateTable() { - if (!TClass::GetClass("St_tpcCorrection")) return 0; -/* - Table: GatingGridB : parameters from Timothy Camarda, 08/09/21 - Old Gating Grid for Run < 18: t0 = 320 ns, setting time = 2.5 us, tau = 2.5 / 4.6 = 543 ns - New Gating Grid for Run 18: t0 = 240 ns, setting time = 1.43 us, tau = 1.43 / 4.6 = 311 ns - Old Gating Grid for Runs 19 - 21 - New Gating Grid for Run > 21: t0 = 240 ns, setting time = 2.0 us, tau = 2.0 / 4.6 = 435 ns. - setting time = time of reaching transperency 99% - tau = setting time/4.6 => exp(-4.6) = 1% - description: Gating Grid transperancy = 0, for t < t0, and 1 - exp(-(t-t0)/tau), for t > t0 -// Fit from 19GeV_2019 data 04/08/2022 G3GFdEdx.root - TF1 *GG = new TF1("GG","TMath::Log(1 - TMath::Exp(-(x-[0])*4.6/[1]))",0,10) - GG->SetParameters(0,2.5) - FitP->Draw("mu:y>>O","i&&j&&abs(x)>40.5","prof") - O->Fit(GG,"er","",0.3,6) - FCN=1174.76 FROM MINOS STATUS=SUCCESSFUL 12 CALLS 109 TOTAL - EDM=1.21976e-08 STRATEGY= 1 ERROR MATRIX ACCURATE - EXT PARAMETER STEP FIRST - NO. NAME VALUE ERROR SIZE DERIVATIVE - 1 p0 1.84098e-01 2.77863e-03 -1.12239e-05 -3.12034e-01 - 2 p1 1.75245e+00 7.48903e-03 7.48903e-03 -3.75507e-04 - 1 p0 2.85771e-01 2.81202e-03 -1.15678e-05 -3.50969e-01 - 2 p1 1.68566e+00 7.58048e-03 7.58048e-03 -4.47722e-04 - - FitP->Draw("mu:y>>I","i&&j&&abs(x)<40.5","prof") -root.exe [72] I->Fit(GG,"er","",0.3,6) - FCN=3798.34 FROM MINOS STATUS=SUCCESSFUL 25 CALLS 119 TOTAL - EDM=1.73946e-09 STRATEGY= 1 ERROR MATRIX ACCURATE - EXT PARAMETER STEP FIRST - NO. NAME VALUE ERROR SIZE DERIVATIVE - 1 p0 1.33585e-01 6.38740e-03 -3.01241e-05 2.44724e-03 - 2 p1 1.29436e+00 1.79433e-02 1.79433e-02 -2.11496e-03 - 1 p0 2.42155e-01 6.17075e-03 -3.49149e-05 2.97485e-03 - 2 p1 1.23716e+00 1.72556e-02 1.72556e-02 -2.32131e-03 -*/ - Int_t nrows = 2; - St_tpcCorrection *tableSet = new St_tpcCorrection("GatingGridB",nrows); - tpcCorrection_st row; - memset(&row,0,tableSet->GetRowSize()); - row.idx = 1; - row.nrows = nrows; - row.npar = 2; - row.a[0] = 2.85771e-01; // t0 - row.a[1] = 1.68566e+00; // settingTime - tableSet->AddAt(&row); // Outer - memset(&row,0,tableSet->GetRowSize()); - row.idx = 2; - row.nrows = nrows; - row.npar = 2; - row.a[0] = 2.42155e-01; // t0 - row.a[1] = 1.23716e+00; // settingTime - tableSet->AddAt(&row); // Outer - return (TDataSet *)tableSet; -} diff --git a/StarDb/Calibrations/tpc/GatingGridB.y2022.C b/StarDb/Calibrations/tpc/GatingGridB.y2022.C new file mode 100644 index 00000000000..5c8d491363b --- /dev/null +++ b/StarDb/Calibrations/tpc/GatingGridB.y2022.C @@ -0,0 +1,25 @@ +#ifndef __CINT__ +#include "tables/St_tpcCorrection_Table.h" +#endif + +TDataSet *CreateTable() { + if (!TClass::GetClass("St_tpcCorrection")) return 0; +/* + Table: GatingGridB : parameters from Timothy Camarda, 08/09/21 + New Gating Grid for Run > 21: t0 = 240 ns, setting time = 2.0 us, tau = 2.0 / 4.6 = 435 ns. + setting time = time of reaching transperency 99% + tau = setting time/4.6 => exp(-4.6) = 1% + description: Gating Grid transperancy = 0, for t < t0, and 1 - exp(-(t-t0)/tau), for t > t0 +*/ + Int_t nrows = 1; + St_tpcCorrection *tableSet = new St_tpcCorrection("GatingGridB",nrows); + tpcCorrection_st row; + memset(&row,0,tableSet->GetRowSize()); + row.idx = 1; + row.nrows = nrows; + row.npar = 2; + row.a[0] = 0.24; // t0 + row.a[1] = 2.00; // settingTime + tableSet->AddAt(&row); // Outer + return (TDataSet *)tableSet; +} diff --git a/StarDb/Calibrations/tpc/starTriggerDelay.C b/StarDb/Calibrations/tpc/starTriggerDelay.C new file mode 120000 index 00000000000..78052916102 --- /dev/null +++ b/StarDb/Calibrations/tpc/starTriggerDelay.C @@ -0,0 +1 @@ +starTriggerDelay.y2019.C \ No newline at end of file diff --git a/StarDb/Geometry/tpc/tpcPadConfig.C b/StarDb/Geometry/tpc/tpcPadConfig.C deleted file mode 100644 index e41376544c0..00000000000 --- a/StarDb/Geometry/tpc/tpcPadConfig.C +++ /dev/null @@ -1,16 +0,0 @@ -#include "tables/St_tpcPadConfig_Table.h" - -TDataSet *CreateTable() { - // ----------------------------------------------------------------- - // db/.const/StarDb/Geometry/tpc/.tpcPadConfig/tpcPadConfig Allocated rows: 1 Used rows: 1 Row size: 1392 bytes - // Table: tpcPadConfig_st[0]--> tpcPadConfig_st[0] - // ==================================================================== - // ------ Test whether this table share library was loaded ------ - if (!TClass::GetClass("St_tpcPadConfig")) return 0; - tpcPadConfig_st row; - St_tpcPadConfig *tableSet = new St_tpcPadConfig("tpcPadConfig",1); - for (Int_t i = 1; i <=24; i++) row.itpc[i-1] = 0; - tableSet->AddAt(&row); - // ----------------- end of code --------------- - return (TDataSet *)tableSet; -} diff --git a/StarDb/Geometry/tpc/tpcPadConfig.y2018.C b/StarDb/Geometry/tpc/tpcPadConfig.y2018.C deleted file mode 100644 index 5bf72970745..00000000000 --- a/StarDb/Geometry/tpc/tpcPadConfig.y2018.C +++ /dev/null @@ -1,19 +0,0 @@ -#include "tables/St_tpcPadConfig_Table.h" - -TDataSet *CreateTable() { - // ----------------------------------------------------------------- - // db/.const/StarDb/Geometry/tpc/.tpcPadConfig/tpcPadConfig Allocated rows: 1 Used rows: 1 Row size: 1392 bytes - // Table: tpcPadConfig_st[0]--> tpcPadConfig_st[0] - // ==================================================================== - // ------ Test whether this table share library was loaded ------ - if (!TClass::GetClass("St_tpcPadConfig")) return 0; - tpcPadConfig_st row; - St_tpcPadConfig *tableSet = new St_tpcPadConfig("tpcPadConfig",1); - Int_t i = 0; - for (i = 1; i <=24; i++) row.itpc[i-1] = 0; - i = 20; - row.itpc[i-1] = 1; - tableSet->AddAt(&row); - // ----------------- end of code --------------- - return (TDataSet *)tableSet; -} diff --git a/StarDb/Geometry/tpc/tpcPadConfig.y2019.C b/StarDb/Geometry/tpc/tpcPadConfig.y2019.C deleted file mode 100644 index a73554331c4..00000000000 --- a/StarDb/Geometry/tpc/tpcPadConfig.y2019.C +++ /dev/null @@ -1,17 +0,0 @@ -#include "tables/St_tpcPadConfig_Table.h" - -TDataSet *CreateTable() { - // ----------------------------------------------------------------- - // db/.const/StarDb/Geometry/tpc/.tpcPadConfig/tpcPadConfig Allocated rows: 1 Used rows: 1 Row size: 1392 bytes - // Table: tpcPadConfig_st[0]--> tpcPadConfig_st[0] - // ==================================================================== - // ------ Test whether this table share library was loaded ------ - if (!TClass::GetClass("St_tpcPadConfig")) return 0; - tpcPadConfig_st row; - St_tpcPadConfig *tableSet = new St_tpcPadConfig("tpcPadConfig",1); - Int_t i = 0; - for (i = 1; i <=24; i++) row.itpc[i-1] = 1; - tableSet->AddAt(&row); - // ----------------- end of code --------------- - return (TDataSet *)tableSet; -} diff --git a/StarDb/dEdxModel/PiDC/pidCorrection.19GeV_2019_P23id.C b/StarDb/dEdxModel/PiDC/pidCorrection.19GeV_2019_P23id.C new file mode 100644 index 00000000000..bc9897535ad --- /dev/null +++ b/StarDb/dEdxModel/PiDC/pidCorrection.19GeV_2019_P23id.C @@ -0,0 +1,50 @@ +TDataSet *CreateTable() { +// ------ Test whether this table share library was loaded ------ + if (!gROOT->GetClass("St_pidCorrection")) return 0; + struct pid_t { // replace a[10] => a0, a1, ..., a9 for cint + Int_t idx; // row index > 0 if it is real + Int_t nrows; // total no. of real rows in the table; For Db interface (where nrows = 50) + Int_t type; // type = 0 polymonical fit, use only [min,max] + Int_t var; // fit variable: 0 => pmomL10, 1 => bgL10, + Int_t particle; // from StEvent/StPidParticleDefinition.h : kUndef = -1, kPidElectron = 0, Proton = 1, Kaon = 2, Pion = 3, Muon = 4, Deuteron = 5, Triton = 6, + // He3 = 7, Alpha = 8, He6 = 9, Li5 = 10, Li6,= 11, Li7 = 12, Be7 = 13, Be9 = 14, Be10 = 15, B11 = 16 + Int_t charge; // +/-1, 0 does not matter + Int_t pull; // != 0 calculated pull correction, == 0 to value + Int_t det; // from StdEdxY2Maker/StTrackPiD.h : kUndef = 0, kI70 = 1, kI70U = 2, kFit = 3, kFitU = 4, kdNdx = 5, kdNdxU = 6, kBTof -7 , kETof = 8, kMtd = 9, kBEmc = 10 + Int_t npar; // npar < 0, X = exp(x) paramterization; abs(npar) >= 100 cut on range [min.max] + Double_t OffSet; // for future use + Double_t min; // fit range + Double_t max; // + Double_t a0; // a[npar] + Double_t a1; // a[npar] + Double_t a2; // a[npar] + Double_t a3; // a[npar] + Double_t a4; // a[npar] + Double_t a5; // a[npar] + Double_t a6; // a[npar] + Double_t a7; // a[npar] + Double_t a8; // a[npar] + Double_t a9; // a[npar] + Char_t comment[32]; + }; + pid_t rows[] = { + //idx,nrows,type,var,part,ch,pull,det,npar,off, min, max, a[10] ,comment + { 0, 0, 0, 0, 3, 0, 0, 3, 4, 0,-0.40, 0.70, 0.0054041, -0.00632595, -0.0149973, -0.00102721, 0,0,0,0,0,0,"PionDEV_dEdx pol3"}, + { 0, 0, 0, 0, 3, 0, 0, 3, 5, 0,-1.10,-0.40, -0.365241, -2.71187, -6.99333, -7.441, -2.76148,0,0,0,0,0,"PionDEV_dEdx pol4"}, + { 0, 0, 0, 0, 1, 0, 0, 3, 4, 0,-0.20, 0.81, 0.0147197, -0.0258435, 0.0812796, -0.140962, 0,0,0,0,0,0,"ProtonDEV_dEdx pol3"}, + { 0, 0, 0, 0, 1, 0, 0, 3, 3, 0,-0.65,-0.20, 0.000859774, -0.110058, 0.0587036, 0,0,0,0,0,0,0,"ProtonDEV_dEdx pol2"}, + { 0, 0, 0, 0, 1, 0, 0, 3, 3, 0,-1.00,-0.65, -1.74183, -5.57528, -4.24693, 0,0,0,0,0,0,0,"ProtonDEV_dEdx pol2"}, + { 0, 0, 0, 0, 0, 0, 0, 3, 4, 0,-0.75, 0.45, 0.121418, 0.0562044, -0.0726388, -0.031316, 0,0,0,0,0,0,"EDEV_dEdx pol3"}, + { 0, 0, 0, 0, 0, 0, 0, 3, 4, 0,-0.85,-0.75, -0.150987, -0.309554, 0.196893, 0.330907, 0,0,0,0,0,0,"EDEV_dEdx pol3"}, + { 0, 0, 0, 0, 0, 0, 0, 3, 4, 0,-1.20,-0.85, -3.9045, -12.4966, -12.8273, -4.23391, 0,0,0,0,0,0,"EDEV_dEdx pol3"} + }; + Int_t nrows = sizeof(rows)/sizeof(pidCorrection_st); + St_pidCorrection *tableSet = new St_pidCorrection("pidCorrection",nrows); + for (Int_t i = 0; i < nrows; i++) { + rows[i].idx = i+1; + rows[i].nrows = nrows; + tableSet->AddAt(&rows[i].idx); + } + // ----------------- end of code --------------- + return (TDataSet *)tableSet; +} diff --git a/StarDb/dEdxModel/PiDC/pidCorrection.C b/StarDb/dEdxModel/PiDC/pidCorrection.C new file mode 100644 index 00000000000..9724719d734 --- /dev/null +++ b/StarDb/dEdxModel/PiDC/pidCorrection.C @@ -0,0 +1,44 @@ +TDataSet *CreateTable() { +// ------ Test whether this table share library was loaded ------ + if (!gROOT->GetClass("St_pidCorrection")) return 0; + struct pid_t { // replace a[10] => a0, a1, ..., a9 for cint + Int_t idx; // row index > 0 if it is real + Int_t nrows; // total no. of real rows in the table; For Db interface (where nrows = 50) + Int_t type; // type = 0 polymonical fit, use only [min,max] + Int_t var; // fit variable: 0 => pmomL10, 1 => bgL10, + Int_t particle; // from StEvent/StPidParticleDefinition.h : kUndef = -1, kPidElectron = 0, Proton = 1, Kaon = 2, Pion = 3, Muon = 4, Deuteron = 5, Triton = 6, + // He3 = 7, Alpha = 8, He6 = 9, Li5 = 10, Li6,= 11, Li7 = 12, Be7 = 13, Be9 = 14, Be10 = 15, B11 = 16 + Int_t charge; // +/-1, 0 does not matter + Int_t pull; // != 0 calculated pull correction, == 0 to value + Int_t det; // from StdEdxY2Maker/StTrackPiD.h : kUndef = 0, kI70 = 1, kI70U = 2, kFit = 3, kFitU = 4, kdNdx = 5, kdNdxU = 6, kBTof -7 , kETof = 8, kMtd = 9, kBEmc = 10 + Int_t npar; // npar < 0, X = exp(x) paramterization; abs(npar) >= 100 cut on range [min.max] + Double_t OffSet; // for future use + Double_t min; // fit range + Double_t max; // + Double_t a0; // a[npar] + Double_t a1; // a[npar] + Double_t a2; // a[npar] + Double_t a3; // a[npar] + Double_t a4; // a[npar] + Double_t a5; // a[npar] + Double_t a6; // a[npar] + Double_t a7; // a[npar] + Double_t a8; // a[npar] + Double_t a9; // a[npar] + Char_t comment[32]; + }; + pid_t rows[] = { + //idx,nrows,type,var,part,ch,pull,det,npar,off,min,max,a[10] + { 0, 0, 0, 0, 0, 0, 0, 3, 0, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,""}, + { 0, 0, 0, 0, 0, 0, 1, 3, 0, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,""} + }; + Int_t nrows = sizeof(rows)/sizeof(pidCorrection_st); + St_pidCorrection *tableSet = new St_pidCorrection("pidCorrection",nrows); + for (Int_t i = 0; i < nrows; i++) { + rows[i].idx = i+1; + rows[i].nrows = nrows; + tableSet->AddAt(&rows[i].idx); + } + // ----------------- end of code --------------- + return (TDataSet *)tableSet; +} diff --git a/StarDb/dEdxModel/PiDC/pidCorrection.r2019.C b/StarDb/dEdxModel/PiDC/pidCorrection.r2019.C new file mode 120000 index 00000000000..d23473efe25 --- /dev/null +++ b/StarDb/dEdxModel/PiDC/pidCorrection.r2019.C @@ -0,0 +1 @@ +pidCorrection.19GeV_2019_P23id.C \ No newline at end of file diff --git a/StarDb/dEdxModel/dNdx_Bichsel.root b/StarDb/dEdxModel/dNdx_Bichsel.root index e197a13cf8a..804184cdfad 100644 Binary files a/StarDb/dEdxModel/dNdx_Bichsel.root and b/StarDb/dEdxModel/dNdx_Bichsel.root differ diff --git a/StarDb/dEdxModel/spline3LndNdxL10.C b/StarDb/dEdxModel/spline3LndNdxL10.C new file mode 100644 index 00000000000..6890be2658f --- /dev/null +++ b/StarDb/dEdxModel/spline3LndNdxL10.C @@ -0,0 +1,13 @@ +TDataSet *CreateTable() { + if (!gROOT->GetClass("St_spline3")) return 0; + Int_t nrows = 1; + St_spline3 *tableSet = new St_spline3("spline3LndNdxL10",nrows); + spline3_st row; + memset(&row,0,tableSet->GetRowSize()); + row.nknots = 14; + Double_t X[14] = {-2.05,-1.5,-1,-0.5,0,0.25,0.5,0.75,1,1.5,2,3,4,5.15}; + Double_t Y[14] = {10.8534,9.17952,7.2945,5.35299,3.83234,3.49396,3.39275,3.41151,3.46255,3.59647,3.67184,3.69744,3.70087,3.70114}; + for (Int_t i = 0; i < 14; i++) {row.Xknots[i] = X[i]; row.Yknots[i] = Y[i];} + tableSet->AddAt(&row); + return (TDataSet *)tableSet; +}