StRoot: StBichsel/dEdxParameterization.cxx Source File

00001 //#define PRINT 1
00002 #define PRINT0 1
00003 //#define PRINT2 1
00004 //#define PRINT3 1
00005 //#define NumDer
00006 #include <assert.h>
00007 #include <Stiostream.h>
00008 #include "dEdxParameterization.h"   
00009 #include "TSystem.h"
00010 #include "TROOT.h"
00011 #include "TDirectory.h"
00012 #include "TF1.h"
00013 #include "TError.h" 
00014 #include "TMath.h"
00015 #define  PrP(A)  cout << "\t" << (#A) << " = \t" << ( A )
00016 
00017 ClassImp(dEdxParameterization)
00018 //________________________________________________________________________________
00019 dEdxParameterization::dEdxParameterization(const Char_t *Tag, Int_t keep3D,
00020                                            const Double_t MostProbableZShift,
00021                                            const Double_t AverageZShift,
00022                                            const Double_t I70Shift,
00023                                            const Double_t I60Shift):
00024   fTag (Tag), fP(0), fA(0), fI70 (0), fI60(0), fD(0),  
00025   fRms (0), fW(0), fPhi(0),
00026   fMostProbableZShift(MostProbableZShift),
00027   fAverageZShift(AverageZShift),
00028   fI70Shift(I70Shift),
00029   fI60Shift(I60Shift),
00030   fbgL10min(-1), fbgL10max(4),
00031   fdxL2min(-0.3), fdxL2max(3),
00032   fzmin(-4), fzmax(6)
00033 {
00034   TDirectory *dir = gDirectory;
00035   const Char_t                                   *rootf = "P10T.root";
00036   if (fTag.Contains("pai" ,TString::kIgnoreCase)) rootf = "PaiT.root";
00037   if (fTag.Contains("p10" ,TString::kIgnoreCase)) rootf = "P10T.root";
00038   if (fTag.Contains("bich",TString::kIgnoreCase)) rootf = "BichselT.root";
00039   static const Char_t *path  = ".:./StarDb/dEdxModel:./StarDb/global/dEdx:./StRoot/StBichsel:$STAR/StarDb/dEdxModel:$STAR/StarDb/global/dEdx:$STAR/StRoot/StBichsel";
00040   Char_t *file = gSystem->Which(path,rootf,kReadPermission);
00041   if (! file) Fatal("dEdxParameterization::GetFile","File %s has not been found in path %s",rootf,path);
00042   else        Warning("dEdxParameterization::GetFile","File %s has been found as %s",rootf,file);
00043   TFile       *pFile = new TFile(file);
00044   assert(pFile);
00045   fP   = (TProfile2D *) pFile->Get("bichP");   assert(fP);   fP->SetDirectory(0);
00046   fA   = (TProfile2D *) pFile->Get("bichA");   assert(fA);   fA->SetDirectory(0);
00047   fI70 = (TProfile2D *) pFile->Get("bichI70"); assert(fI70); fI70->SetDirectory(0);
00048   fI60 = (TProfile2D *) pFile->Get("bichI60"); assert(fI60); fI60->SetDirectory(0);
00049   fD   = (TProfile2D *) pFile->Get("bichD");   assert(fD);   fD->SetDirectory(0);
00050   fRms = (TProfile2D *) pFile->Get("bichRms"); assert(fRms); fRms->SetDirectory(0);
00051   fW   = (TProfile2D *) pFile->Get("bichW");   assert(fW);   fW->SetDirectory(0);
00052   fPhi = (TH3D       *) pFile->Get("bichPhi"); assert(fPhi); fPhi->SetDirectory(0);
00053   delete pFile;  
00054   fbgL10min = fPhi->GetXaxis()->GetBinCenter(1) + 1e-7;
00055   fbgL10max = fPhi->GetXaxis()->GetBinCenter(fPhi->GetXaxis()->GetNbins()) - 1e-7;
00056   fdxL2min  = fPhi->GetYaxis()->GetBinCenter(1) + 1e-7;
00057   fdxL2max  = fPhi->GetYaxis()->GetBinCenter(fPhi->GetYaxis()->GetNbins()) - 1e-7;
00058   fzmin  = fPhi->GetZaxis()->GetBinCenter(1) + 1e-7;
00059   fzmax  = fPhi->GetZaxis()->GetBinCenter(fPhi->GetZaxis()->GetNbins()) - 1e-7;
00060   if (dir) dir->cd();
00061   for (Int_t i = 0; i<3; i++) {
00062     if (i == 0) fAXYZ[i] = fPhi->GetXaxis(); 
00063     if (i == 1) fAXYZ[i] = fPhi->GetYaxis(); 
00064     if (i == 2) fAXYZ[i] = fPhi->GetZaxis(); 
00065     fnBins[i] = fAXYZ[i]->GetNbins(); 
00066     fbinW[i]  = fAXYZ[i]->GetBinWidth(1);
00067 #ifdef PRINT0
00068     PrP(i); PrP(fnBins[i]); PrP(fbinW[i]); cout << endl;
00069     assert(fnBins[i] !=1);
00070 #endif
00071   }
00072   //  if (! keep3D) SafeDelete(fPhi);
00073   // set normalization factor to 2.3976 keV/cm at beta*gamma = 4;
00074   static const Double_t dEdxMIP = 2.39761562607903311; // [keV/cm]
00075   static const Double_t MIPBetaGamma10 = TMath::Log10(4.);
00076   //  fMostProbableZShift = TMath::Log(dEdxMIP) - Interpolation(fP,MIPBetaGamma10,1,0);
00077   //  fAverageZShift      = TMath::Log(dEdxMIP) - Interpolation(fA,MIPBetaGamma10,1,0);
00078   fI70Shift           *= dEdxMIP/GetI70(MIPBetaGamma10,1);
00079   fI60Shift           *= dEdxMIP/GetI60(MIPBetaGamma10,1);
00080   fMostProbableZShift  = TMath::Log(fI70Shift);
00081   fAverageZShift       = fMostProbableZShift;
00082 }
00083 //________________________________________________________________________________
00084 dEdxParameterization::~dEdxParameterization() { 
00085   SafeDelete(fP);  
00086   SafeDelete(fA);  
00087   SafeDelete(fI70);
00088   SafeDelete(fI60);
00089   SafeDelete(fD);  
00090   SafeDelete(fRms);
00091   SafeDelete(fW);
00092   SafeDelete(fPhi);
00093 }    
00094 #if ROOT_VERSION_CODE < 334336   /* 5.26.0 */
00095 //________________________________________________________________________________
00096 Double_t    dEdxParameterization::Interpolation(Int_t Narg, TH1 *hist, Double_t *XYZ, Int_t kase) {
00097   assert(hist);
00098   Int_t Ndim = hist->GetDimension();
00099   assert (Ndim>=1 && Ndim <= 3);
00100   assert (Ndim == Narg);
00101 #if defined(PRINT) || defined(PRINT2)
00102   cout << "Interpolation:";
00103   PrP(XYZ[0]); PrP(XYZ[1]); PrP(XYZ[2]); 
00104   cout << endl;
00105 #endif
00106   Double_t Value = 0;
00107   Int_t iXYZ[3], ixyz[3];
00108   Double_t dXYZ[3], pXYZ[3];
00109   Int_t i;
00110   for (i = 0; i< Ndim; i++) {
00111     iXYZ[i] = fAXYZ[i]->FindBin(XYZ[i]); 
00112     //    if (iXYZ[i] <= 0 || iXYZ[i] >= fnBins[i]) return Value; // outside of the range
00113     if (iXYZ[i] < 2) iXYZ[i] = 2; 
00114     if (iXYZ[i] >= fnBins[i])  iXYZ[i] =  fnBins[i] - 1; 
00115     dXYZ[i] = (XYZ[i] - fAXYZ[i]->GetBinCenter(iXYZ[i]))/fbinW[i];
00116     if (dXYZ[i] >=0) ixyz[i] = iXYZ[i] + 1;
00117     else            {ixyz[i] = iXYZ[i] - 1; dXYZ[i] = - dXYZ[i];}
00118     pXYZ[i] = 1. - dXYZ[i];
00119     if (i == kase - 1) {// derivatives wrt X
00120       dXYZ[i] = 1./fbinW[i];
00121       if (ixyz[i] < iXYZ[i]) dXYZ[i] = - dXYZ[i];
00122       pXYZ[i] =    -dXYZ[i]; 
00123     }
00124   } 
00125   
00126 #ifdef PRINT
00127   for (i = 0; i<Ndim; i++) {PrP(i); PrP(iXYZ[i]);} cout << endl;
00128   for (i = 0; i<Ndim; i++) {PrP(i); PrP(dXYZ[i]);} cout << endl;
00129   for (i = 0; i<Ndim; i++) {PrP(i); PrP(pXYZ[i]);} cout << endl;
00130   
00131   PrP(iXYZ[0]);PrP(iXYZ[1]);PrP(iXYZ[2]); PrP(hist->GetBinContent(iXYZ[0],iXYZ[1],iXYZ[2])); cout << endl; 
00132   PrP(ixyz[0]);PrP(iXYZ[1]);PrP(iXYZ[2]); PrP(hist->GetBinContent(ixyz[0],iXYZ[1],iXYZ[2])); cout << endl; 
00133   PrP(iXYZ[0]);PrP(ixyz[1]);PrP(iXYZ[2]); PrP(hist->GetBinContent(iXYZ[0],ixyz[1],iXYZ[2])); cout << endl; 
00134   PrP(ixyz[0]);PrP(ixyz[1]);PrP(iXYZ[2]); PrP(hist->GetBinContent(ixyz[0],ixyz[1],iXYZ[2])); cout << endl;
00135   PrP(iXYZ[0]);PrP(iXYZ[1]);PrP(ixyz[2]); PrP(hist->GetBinContent(iXYZ[0],iXYZ[1],ixyz[2])); cout << endl; 
00136   PrP(ixyz[0]);PrP(iXYZ[1]);PrP(ixyz[2]); PrP(hist->GetBinContent(ixyz[0],iXYZ[1],ixyz[2])); cout << endl; 
00137   PrP(iXYZ[0]);PrP(ixyz[1]);PrP(ixyz[2]); PrP(hist->GetBinContent(iXYZ[0],ixyz[1],ixyz[2])); cout << endl; 
00138   PrP(ixyz[0]);PrP(ixyz[1]);PrP(ixyz[2]); PrP(hist->GetBinContent(ixyz[0],ixyz[1],ixyz[2])); cout << endl;
00139 #endif  
00140   switch (Ndim) {
00141   case 3:  Value = 
00142              pXYZ[0]*pXYZ[1]*pXYZ[2]*hist->GetBinContent(iXYZ[0],iXYZ[1],iXYZ[2]) + 
00143              dXYZ[0]*pXYZ[1]*pXYZ[2]*hist->GetBinContent(ixyz[0],iXYZ[1],iXYZ[2]) + 
00144              pXYZ[0]*dXYZ[1]*pXYZ[2]*hist->GetBinContent(iXYZ[0],ixyz[1],iXYZ[2]) + 
00145              dXYZ[0]*dXYZ[1]*pXYZ[2]*hist->GetBinContent(ixyz[0],ixyz[1],iXYZ[2]) +
00146              pXYZ[0]*pXYZ[1]*dXYZ[2]*hist->GetBinContent(iXYZ[0],iXYZ[1],ixyz[2]) + 
00147              dXYZ[0]*pXYZ[1]*dXYZ[2]*hist->GetBinContent(ixyz[0],iXYZ[1],ixyz[2]) + 
00148              pXYZ[0]*dXYZ[1]*dXYZ[2]*hist->GetBinContent(iXYZ[0],ixyz[1],ixyz[2]) + 
00149              dXYZ[0]*dXYZ[1]*dXYZ[2]*hist->GetBinContent(ixyz[0],ixyz[1],ixyz[2]);
00150   break;
00151   case 2: Value = 
00152             pXYZ[0]*pXYZ[1]*hist->GetBinContent(iXYZ[0],iXYZ[1]) + 
00153             dXYZ[0]*pXYZ[1]*hist->GetBinContent(ixyz[0],iXYZ[1]) + 
00154             pXYZ[0]*dXYZ[1]*hist->GetBinContent(iXYZ[0],ixyz[1]) + 
00155             dXYZ[0]*dXYZ[1]*hist->GetBinContent(ixyz[0],ixyz[1]);
00156   break;
00157   case 1: Value = 
00158             pXYZ[0]*hist->GetBinContent(iXYZ[0]) + dXYZ[0]*hist->GetBinContent(ixyz[0]);
00159   break;
00160   default:
00161     assert(0);
00162   }
00163 #if defined(PRINT) || defined(PRINT2)
00164   PrP(Value); cout << endl;
00165 #endif
00166   return Value;
00167 }
00168 //________________________________________________________________________________
00169 Double_t    dEdxParameterization::Interpolation(TH3 *hist, Double_t X, Double_t Y, Double_t Z, Int_t kase) {
00170   assert(hist);
00171 #if defined(PRINT) || defined(PRINT2)
00172   cout << "Interpolation:";
00173   PrP(X); PrP(Y); PrP(Z); 
00174   cout << endl;
00175 #endif
00176   Double_t XYZ[3];
00177   XYZ[0] = X;
00178   XYZ[1] = Y;
00179   XYZ[2] = Z;
00180   return Interpolation(3, hist, XYZ, kase);
00181 }
00182 #endif /* ROOT_VERSION_CODE < ROOT_VERSION(5,26,0)  */
00183 //________________________________________________________________________________
00184 void dEdxParameterization::Print() {
00185   PrP(fTag); cout << endl;
00186   PrP(fP); if (fP) PrP(fP->GetTitle()); cout << endl;
00187   PrP(fA); if (fA) PrP(fA->GetTitle()); cout << endl;
00188   PrP(fI70); if (fI70) PrP(fI70->GetTitle()); cout << endl;
00189   PrP(fI60); if (fI60) PrP(fI60->GetTitle()); cout << endl;
00190   PrP(fD); if (fD) PrP(fD->GetTitle()); cout << endl;
00191   PrP(fRms); if (fRms) PrP(fRms->GetTitle()); cout << endl;
00192   PrP(fW); if (fW) PrP(fW->GetTitle()); cout << endl;
00193   PrP(fPhi); if (fPhi) PrP(fPhi->GetTitle()); cout << endl;
00194   PrP(fMostProbableZShift); cout << endl;
00195   PrP(fAverageZShift); cout << endl;
00196   PrP(fI70Shift); cout << endl;
00197   PrP(fI60Shift); cout << endl;
00198 } 
00199 //________________________________________________________________________________
00200 Double_t dEdxParameterization::MostProbableZCorrection(Double_t log10bg) {
00201   static const Double_t pars[2] = {-3.68846e-03, 4.72944e+00}; // FitHzAllHist012P05id  FitH + Prof 050905
00202   return pars[0]*TMath::Exp(-pars[1]*log10bg);
00203 }//________________________________________________________________________________
00204 Double_t dEdxParameterization::I70Correction(Double_t log10bg) {
00205   static const Double_t pars[2] = {-1.65714e-02, 3.27271e+00}; //  FitH70AllHist012P05id FitH + Prof 050905
00206   return TMath::Exp(pars[0]*TMath::Exp(-pars[1]*log10bg));
00207 }