StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
dProjector.cxx
1 /***************************************************************************
2  *
3  * Author: Dominik Flierl, flierl@bnl.gov
4  ***************************************************************************
5  *
6  * Description: part of STAR HBT Framework: StHbtMaker package
7  * this is an object which produces projections of 3d histograms
8  *
9  **************************************************************************/
10 
11 #include "dProjector.h"
12 #include "TVector3.h"
13 
14 #ifdef __ROOT__
15 ClassImp(dProjector)
16 #endif
17 
18 #define __ABS__(x) (x>0) ? x : -x
19 
20 
21 //___________________
22 dProjector::dProjector(dFitter3d* dFitter)
23 {
24  // get the 3d histo
25  mRatio = dFitter->Ratio() ;
26 
27  // if the ratio histogram does not exist -> generate it
28  if (!mRatio)
29  {
30  mRatio = (TH3D*) (dFitter->Numerator())->Clone() ;
31  mRatio->SetName("ratio") ;
32  mRatio->Scale(1.0/dFitter->Norm()) ;
33  mRatio->Divide(dFitter->Denominator()) ;
34  }
35 
36  // get minuit to extract the fit result
37  TMinuit* tMinuit = dFitter->getMinuit() ;
38  // fill fit result into array
39  mFitParameters = new double[5] ;
40  cout << " used fit parameters : " << "\t" ;
41  for( int index = 0 ; index < 5 ; index++ )
42  {
43  double value = 0 ;
44  double error = 0 ;
45  tMinuit->GetParameter(index, value, error) ;
46  mFitParameters[index] = value ;
47  cout << value << "\t" ;
48  }
49  cout << endl ;
50 } ;
51 //___________________
52 dProjector::~dProjector()
53 {
54  delete [] mFitParameters ;
55 } ;
56 //___________________
57 TH1D* dProjector::get1dProjection(TString axis, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax)
58 {
59  // get ranges of the ratio 3d histogram
60  int xminBin = mRatio->GetXaxis()->FindFixBin(xmin) ;
61  int xmaxBin = mRatio->GetXaxis()->FindFixBin(xmax) ;
62  int yminBin = mRatio->GetYaxis()->FindFixBin(ymin) ;
63  int ymaxBin = mRatio->GetYaxis()->FindFixBin(ymax) ;
64  int zminBin = mRatio->GetZaxis()->FindFixBin(zmin) ;
65  int zmaxBin = mRatio->GetZaxis()->FindFixBin(zmax) ;
66 
67  // get dimensions of the projection histos
68  int proNbins = 0 ;
69  double proMin = 0.0 ;
70  double proMax = 0.0 ;
71  char* axistitle = new char[20] ;
72  if (axis == "x") { proNbins = __ABS__(xmaxBin-xminBin)+1 ; proMin = xmin ; proMax = xmax ; strcpy(axistitle,mRatio->GetXaxis()->GetTitle()) ; }
73  else if (axis == "y") { proNbins = __ABS__(ymaxBin-yminBin)+1 ; proMin = ymin ; proMax = ymax ; strcpy(axistitle,mRatio->GetYaxis()->GetTitle()) ; }
74  else if (axis == "z") { proNbins = __ABS__(zmaxBin-zminBin)+1 ; proMin = zmin ; proMax = zmax ; strcpy(axistitle,mRatio->GetZaxis()->GetTitle()) ; }
75  else { cout << "Wrong axis in projector, this should not happen ! "<< endl ; return 0; } ;
76 
77  // create 1d histo
78  int ran = (int) (xmin*ymin*zmin*10000) ; // to distinguish between objects give them a special random id ...
79  m1dpro = new TH1D("my1dpro","my1dpro",proNbins,proMin,proMax) ;
80  m1dpro->SetXTitle(axistitle) ;
81  m1dpro->SetTitle(axistitle) ;
82  char restitle[20] ; sprintf(restitle,"fit%s%d\n",axis.Data(),ran) ;
83  m1dpro->SetName(restitle) ;
84  m1dfit = new TH1D("my1dfit","my1dfit",proNbins,proMin,proMax) ;
85  char fittitle[20] ; sprintf(fittitle,"fit%s%d\n",axis.Data(),ran) ;
86  m1dfit->SetName(fittitle);
87  m1dfit->SetLineColor(2) ;
88  m1dfit->SetLineStyle(2) ;
89  m1dfit->SetLineWidth(5) ;
90  m1dnor = new TH1D("my1dnor","my1dnor",proNbins,proMin,proMax) ;
91  char normtitle[20] ; sprintf(normtitle,"norm%s%d\n",axis.Data(),ran) ;
92  m1dnor->SetName(normtitle);
93 
94  // fill histo
95  for(int indexX = xminBin ; indexX <= xmaxBin ; indexX++)
96  {
97  for(int indexY = yminBin ; indexY <= ymaxBin ; indexY++)
98  {
99  for(int indexZ = zminBin ; indexZ <= zmaxBin ; indexZ++)
100  {
101  // get the content of the 3d cell
102  double cellcontent = mRatio->GetBinContent(indexX,indexY,indexZ) ;
103  // get the content of the 3d cell using the fit function
104  TVector3 position(mRatio->GetXaxis()->GetBinCenter(indexX),
105  mRatio->GetYaxis()->GetBinCenter(indexY),
106  mRatio->GetZaxis()->GetBinCenter(indexZ));
107  double cellContentFromFit = mFitter->ykpCorrelationFunction(position,mFitParameters) ;
108 
109 
110  // fill the 1d histos
111  int proBin = - 100 ;
112  if (axis == "x") { proBin = m1dpro->GetXaxis()->FindFixBin(mRatio->GetXaxis()->GetBinCenter(indexX)); }
113  else if (axis == "y") { proBin = m1dpro->GetXaxis()->FindFixBin(mRatio->GetYaxis()->GetBinCenter(indexY)); }
114  else if (axis == "z") { proBin = m1dpro->GetXaxis()->FindFixBin(mRatio->GetZaxis()->GetBinCenter(indexZ)); }
115 
116  if(cellcontent>0.5 && cellcontent <1.5)
117  {
118  m1dpro->AddBinContent(proBin,cellcontent) ;
119  m1dfit->AddBinContent(proBin,cellContentFromFit) ;
120 
121  if (axis == "z" && proBin<10) cout << proBin << " " ;
122  // count the number of fillings to normalize the projection to 1
123  m1dnor->AddBinContent(proBin,1.0) ;
124  }
125  }
126  }
127  }
128  // normalize fit and projection histos
129  m1dpro->Divide(m1dnor) ;
130  m1dfit->Divide(m1dnor) ;
131 
132  // return projection
133  return m1dpro ;
134 }
135 //___________________
136 TH2D* dProjector::get2dProjection(TString axis, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax)
137 {
138  // get ranges of the ratio 3d histogram
139  int xminBin = mRatio->GetXaxis()->FindFixBin(xmin) ;
140  int xmaxBin = mRatio->GetXaxis()->FindFixBin(xmax) ;
141  int yminBin = mRatio->GetYaxis()->FindFixBin(ymin) ;
142  int ymaxBin = mRatio->GetYaxis()->FindFixBin(ymax) ;
143  int zminBin = mRatio->GetZaxis()->FindFixBin(zmin) ;
144  int zmaxBin = mRatio->GetZaxis()->FindFixBin(zmax) ;
145 
146  // get dimensions of the projection histos
147  int proNbinsX = 0 ;
148  double proMinX = 0.0 ;
149  double proMaxX = 0.0 ;
150  char* axistitleX = new char[20] ;
151  int proNbinsY = 0 ;
152  double proMinY = 0.0 ;
153  double proMaxY = 0.0 ;
154  char* axistitleY = new char[20] ;
155 
156  if (axis == "x")
157  {
158  proNbinsX = __ABS__(ymaxBin-yminBin)+1 ; proMinX = ymin ; proMaxX = ymax ; strcpy(axistitleX,mRatio->GetYaxis()->GetTitle()) ;
159  proNbinsY = __ABS__(zmaxBin-zminBin)+1 ; proMinY = zmin ; proMaxY = zmax ; strcpy(axistitleY,mRatio->GetZaxis()->GetTitle()) ;
160  }
161  else if (axis == "y")
162  {
163  proNbinsX = __ABS__(zmaxBin-zminBin)+1 ; proMinX = zmin ; proMaxX = zmax ; strcpy(axistitleX,mRatio->GetZaxis()->GetTitle()) ;
164  proNbinsY = __ABS__(xmaxBin-xminBin)+1 ; proMinY = xmin ; proMaxY = xmax ; strcpy(axistitleY,mRatio->GetXaxis()->GetTitle()) ;
165  }
166  else if (axis == "z")
167  {
168  proNbinsX = __ABS__(xmaxBin-xminBin)+1 ; proMinX = xmin ; proMaxX = xmax ; strcpy(axistitleX,mRatio->GetXaxis()->GetTitle()) ;
169  proNbinsY = __ABS__(ymaxBin-yminBin)+1 ; proMinY = ymin ; proMaxY = ymax ; strcpy(axistitleY,mRatio->GetYaxis()->GetTitle()) ;
170  }
171 
172  else { cout << "Wrong axis in projector, this should not happen ! "<< endl ; return 0; } ;
173 
174  // create 2d histo
175  int ran = (int) (xmin*ymin*zmin*10000) ; // to distinguish between objects give them a special random id ...
176  // projection
177  m2dpro = new TH2D("my2dpro","my2dpro",proNbinsX,proMinX,proMaxX,proNbinsY,proMinY,proMaxY) ;
178  char resName[20] ; sprintf(resName,"res%s%d\n",axis.Data(),ran) ; m2dpro->SetName(resName) ;
179  m2dpro->SetXTitle(axistitleX) ;
180  m2dpro->SetYTitle(axistitleY) ;
181  char restitle[20] ; sprintf(restitle,"%s\n",axis.Data()) ;m2dpro->SetTitle(restitle) ;
182  // fit
183  m2dfit = new TH2D("my2dfit","my2dfit",proNbinsX,proMinX,proMaxX,proNbinsY,proMinY,proMaxY) ;
184  char fittitle[20] ; sprintf(fittitle,"fit%s%d\n",axis.Data(),ran) ; m2dfit->SetName(fittitle) ;
185  // norm
186  m2dnor = new TH2D("my2dnor","my2dnor",proNbinsX,proMinX,proMaxX,proNbinsY,proMinY,proMaxY) ;
187  char normtitle[20] ; sprintf(normtitle,"norm%s%d\n",axis.Data(),ran) ; m2dnor->SetName(normtitle) ;
188 
189  // fill histo
190  for(int indexX = xminBin ; indexX <= xmaxBin ; indexX++)
191  {
192  for(int indexY = yminBin ; indexY <= ymaxBin ; indexY++)
193  {
194  for(int indexZ = zminBin ; indexZ <= zmaxBin ; indexZ++)
195  {
196  // get the content of the 3d cell
197  double cellcontent = mRatio->GetBinContent(indexX,indexY,indexZ) ;
198  // get the content of the 3d cell using the fit function
199  TVector3 position(mRatio->GetXaxis()->GetBinCenter(indexX),
200  mRatio->GetYaxis()->GetBinCenter(indexY),
201  mRatio->GetZaxis()->GetBinCenter(indexZ));
202  double cellContentFromFit = mFitter->ykpCorrelationFunction(position,mFitParameters) ;
203 
204 
205  // fill the 1d histos
206  double proX = 100000.0 ;
207  double proY = 100000.0 ;
208  if (axis == "x")
209  {
210  proX = mRatio->GetYaxis()->GetBinCenter(indexY) ;
211  proY = mRatio->GetZaxis()->GetBinCenter(indexZ) ;
212  }
213  else if (axis == "y")
214  {
215  proX = mRatio->GetZaxis()->GetBinCenter(indexZ) ;
216  proY = mRatio->GetXaxis()->GetBinCenter(indexX) ;
217  }
218  else if (axis == "z")
219  {
220  proX = mRatio->GetXaxis()->GetBinCenter(indexX) ;
221  proY = mRatio->GetYaxis()->GetBinCenter(indexY) ;
222  }
223 
224  if(cellcontent>0.5 && cellcontent <1.5)
225  {
226  // now fill 2d histos
227  m2dpro->Fill(proX,proY,cellcontent) ;
228  m2dfit->Fill(proX,proY,cellContentFromFit) ;
229  // count the number of fillings to normalize the projection to 1
230  m2dnor->Fill(proX,proY,1.0) ;
231  }
232  }
233  }
234  }
235  // normalize fit and projection histos
236  m2dpro->Divide(m2dnor) ;
237  m2dfit->Divide(m2dnor) ;
238 
239  // return projection
240  return m2dpro ;
241 
242 }
Generated by   doxygen 1.8.5