StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
dEdxParameterization.cxx
1 //#define PRINT 1
2 #define PRINT0 1
3 //#define PRINT2 1
4 //#define PRINT3 1
5 //#define NumDer
6 #include <assert.h>
7 #include <Stiostream.h>
8 #include "dEdxParameterization.h"
9 #include "TSystem.h"
10 #include "TROOT.h"
11 #include "TDirectory.h"
12 #include "TF1.h"
13 #include "TError.h"
14 #include "TMath.h"
15 
16 #define PrP(A) cout << "\t" << (#A) << " = \t" << ( A )
17 
18 ClassImp(dEdxParameterization)
19 //________________________________________________________________________________
20 dEdxParameterization::dEdxParameterization(const Char_t *Tag, Int_t keep3D,
21  const Double_t MostProbableZShift,
22  const Double_t AverageZShift,
23  const Double_t I70Shift,
24  const Double_t I60Shift):
25  fTag (Tag), fP(0), fA(0), fI70 (0), fI60(0), fD(0),
26  fRms (0), fW(0), fPhi(0),
27  fMostProbableZShift(MostProbableZShift),
28  fAverageZShift(AverageZShift),
29  fI70Shift(I70Shift),
30  fI60Shift(I60Shift),
31  fbgL10min(-1), fbgL10max(4),
32  fdxL2min(-0.3), fdxL2max(3),
33  fzmin(-4), fzmax(6)
34 {
35  memset (fTrs, 0, sizeof(fTrs));
36  TDirectory *dir = gDirectory;
37  const Char_t *rootf = "P10T.root";
38 #if 0
39  if (fTag.Contains("pai" ,TString::kIgnoreCase)) rootf = "PaiT.root";
40  if (fTag.Contains("p10" ,TString::kIgnoreCase)) rootf = "P10T.root";
41  if (fTag.Contains("bich",TString::kIgnoreCase)) rootf = "BichselT.root";
42 #endif
43  static const Char_t *path = ".:./StarDb/dEdxModel:./StarDb/global/dEdx:./StRoot/StBichsel:$STAR/StarDb/dEdxModel:$STAR/StarDb/global/dEdx:$STAR/StRoot/StBichsel";
44  Char_t *file = gSystem->Which(path,rootf,kReadPermission);
45  if (! file) Fatal("dEdxParameterization::GetFile","File %s has not been found in path %s",rootf,path);
46  else Warning("dEdxParameterization::GetFile","File %s has been found as %s",rootf,file);
47  TFile *pFile = new TFile(file);
48  delete [] file;
49  assert(pFile);
50  fP = (TProfile2D *) pFile->Get("bichP"); assert(fP); fP->SetDirectory(0);
51  fA = (TProfile2D *) pFile->Get("bichA"); assert(fA); fA->SetDirectory(0);
52  fI70 = (TProfile2D *) pFile->Get("bichI70"); assert(fI70); fI70->SetDirectory(0);
53  fI60 = (TProfile2D *) pFile->Get("bichI60"); assert(fI60); fI60->SetDirectory(0);
54  fD = (TProfile2D *) pFile->Get("bichD"); assert(fD); fD->SetDirectory(0);
55  fRms = (TProfile2D *) pFile->Get("bichRms"); assert(fRms); fRms->SetDirectory(0);
56  fW = (TProfile2D *) pFile->Get("bichW"); assert(fW); fW->SetDirectory(0);
57  fPhi = (TH3D *) pFile->Get("bichPhi"); assert(fPhi); fPhi->SetDirectory(0);
58  fbgL10min = fPhi->GetXaxis()->GetBinCenter(1) + 1e-7;
59  fbgL10max = fPhi->GetXaxis()->GetBinCenter(fPhi->GetXaxis()->GetNbins()) - 1e-7;
60  fdxL2min = fPhi->GetYaxis()->GetBinCenter(1) + 1e-7;
61  fdxL2max = fPhi->GetYaxis()->GetBinCenter(fPhi->GetYaxis()->GetNbins()) - 1e-7;
62  fzmin = fPhi->GetZaxis()->GetBinCenter(1) + 1e-7;
63  fzmax = fPhi->GetZaxis()->GetBinCenter(fPhi->GetZaxis()->GetNbins()) - 1e-7;
64  if (dir) dir->cd();
65  for (Int_t i = 0; i<3; i++) {
66  if (i == 0) fAXYZ[i] = fPhi->GetXaxis();
67  if (i == 1) fAXYZ[i] = fPhi->GetYaxis();
68  if (i == 2) fAXYZ[i] = fPhi->GetZaxis();
69  fnBins[i] = fAXYZ[i]->GetNbins();
70  fbinW[i] = fAXYZ[i]->GetBinWidth(1);
71 #ifdef PRINT0
72  PrP(i); PrP(fnBins[i]); PrP(fbinW[i]); cout << endl;
73  assert(fnBins[i] !=1);
74 #endif
75  }
76  // if (! keep3D) SafeDelete(fPhi);
77  // set normalization factor to 2.3976 keV/cm at beta*gamma = 4;
78  static const Double_t dEdxMIP = 2.39761562607903311; // [keV/cm]
79  static const Double_t MIPBetaGamma10 = TMath::Log10(4.);
80  // fMostProbableZShift = TMath::Log(dEdxMIP) - Interpolation(fP,MIPBetaGamma10,1,0);
81  // fAverageZShift = TMath::Log(dEdxMIP) - Interpolation(fA,MIPBetaGamma10,1,0);
82  fI70Shift *= dEdxMIP/GetI70(MIPBetaGamma10,1);
83  fI60Shift *= dEdxMIP/GetI60(MIPBetaGamma10,1);
84  fMostProbableZShift = TMath::Log(fI70Shift);
85  fAverageZShift = fMostProbableZShift;
86  const Char_t *Names[KPidParticles+1] = {"e","proton","kaon","pi","mu","deuteron","triton","He3","alpha","all"};
87  for (Int_t i = 0; i <= KPidParticles; i++) {
88  TString name(Names[i]);
89  const Char_t *type[6] = {"70p","70","70S","zp","z","zS"};
90  for (Int_t j = 0; j < 6; j++) {
91  fTrs[i][j] = (TH1D *) pFile->Get(name + type[j]);
92  if (fTrs[i][j]) fTrs[i][j]->SetDirectory(0);
93  }
94  }
95  delete pFile;
96 }
97 //________________________________________________________________________________
98 dEdxParameterization::~dEdxParameterization() {
99  SafeDelete(fP);
100  SafeDelete(fA);
101  SafeDelete(fI70);
102  SafeDelete(fI60);
103  SafeDelete(fD);
104  SafeDelete(fRms);
105  SafeDelete(fW);
106  SafeDelete(fPhi);
107  for (Int_t i = 0; i <= KPidParticles; i++)
108  for (Int_t j = 0; j < 6; j++) {SafeDelete(fTrs[i][j]);}
109 }
110 //________________________________________________________________________________
111 void dEdxParameterization::Print() {
112  PrP(fTag); cout << endl;
113  PrP(fP); if (fP) PrP(fP->GetTitle()); cout << endl;
114  PrP(fA); if (fA) PrP(fA->GetTitle()); cout << endl;
115  PrP(fI70); if (fI70) PrP(fI70->GetTitle()); cout << endl;
116  PrP(fI60); if (fI60) PrP(fI60->GetTitle()); cout << endl;
117  PrP(fD); if (fD) PrP(fD->GetTitle()); cout << endl;
118  PrP(fRms); if (fRms) PrP(fRms->GetTitle()); cout << endl;
119  PrP(fW); if (fW) PrP(fW->GetTitle()); cout << endl;
120  PrP(fPhi); if (fPhi) PrP(fPhi->GetTitle()); cout << endl;
121  PrP(fMostProbableZShift); cout << endl;
122  PrP(fAverageZShift); cout << endl;
123  PrP(fI70Shift); cout << endl;
124  PrP(fI60Shift); cout << endl;
125 }
126 //________________________________________________________________________________
127 Double_t dEdxParameterization::MostProbableZCorrection(Double_t log10bg) {
128  static const Double_t pars[2] = {-3.68846e-03, 4.72944e+00}; // FitHzAllHist012P05id FitH + Prof 050905
129  return pars[0]*TMath::Exp(-pars[1]*log10bg);
130 }//________________________________________________________________________________
131 Double_t dEdxParameterization::I70Correction(Double_t log10bg) {
132  static const Double_t pars[2] = {-1.65714e-02, 3.27271e+00}; // FitH70AllHist012P05id FitH + Prof 050905
133  return TMath::Exp(pars[0]*TMath::Exp(-pars[1]*log10bg));
134 }
135 //________________________________________________________________________________
136 Double_t dEdxParameterization::Get(const TH1D *hist, Double_t log10bg) const {
137  static TH1D *hsave = 0;
138  static Double_t xmin = -100, xmax = 100;
139  if (hist != hsave) {
140  hsave = (TH1D *) hist;
141  TAxis *x = hsave->GetXaxis();
142  Int_t f = x->GetFirst();
143  Int_t l = x->GetLast();
144 
145  xmin = x->GetBinUpEdge(f);
146  xmax = x->GetBinLowEdge(l);
147  }
148  if (log10bg < xmin) log10bg = xmin;
149  if (log10bg > xmax) log10bg = xmax;
150  return hsave->Interpolate(log10bg);
151 }
152 // $Id: dEdxParameterization.cxx,v 1.19 2016/06/10 19:56:11 fisyak Exp $
153 // $Log: dEdxParameterization.cxx,v $
154 // Revision 1.19 2016/06/10 19:56:11 fisyak
155 // Fix covertry warning
156 //
157 // Revision 1.18 2015/12/24 00:16:26 fisyak
158 // Add TpcRS model and macros
159 //
Generated by   doxygen 1.8.5