StRoot: StHbtMaker/Fit/dProjector.cxx Source File

00001 /***************************************************************************
00002  *
00003  * Author: Dominik Flierl, flierl@bnl.gov
00004  ***************************************************************************
00005  *
00006  * Description: part of STAR HBT Framework: StHbtMaker package
00007  *   this is an object which produces projections of 3d histograms
00008  *  
00009  **************************************************************************/
00010 
00011 #include "dProjector.h"
00012 #include "TVector3.h"
00013 
00014 #ifdef __ROOT__
00015 ClassImp(dProjector)
00016 #endif
00017 
00018 #define __ABS__(x) (x>0) ? x : -x
00019 
00020   
00021 //___________________
00022 dProjector::dProjector(dFitter3d* dFitter) 
00023 { 
00024   // get the 3d histo
00025   mRatio = dFitter->Ratio() ;
00026   
00027   // if the ratio histogram does not exist -> generate it
00028   if (!mRatio)
00029     {
00030       mRatio = (TH3D*) (dFitter->Numerator())->Clone() ;
00031       mRatio->SetName("ratio") ;
00032       mRatio->Scale(1.0/dFitter->Norm()) ;
00033       mRatio->Divide(dFitter->Denominator()) ;
00034     }
00035 
00036   // get minuit to extract the fit result
00037   TMinuit* tMinuit = dFitter->getMinuit() ;
00038   // fill fit result into array
00039   mFitParameters = new double[5] ;
00040   cout << " used fit parameters : " << "\t" ;
00041   for( int index = 0 ; index < 5 ; index++ )
00042     {
00043       double value = 0 ;
00044       double error = 0 ;
00045       tMinuit->GetParameter(index, value, error) ;
00046       mFitParameters[index] = value ;
00047       cout << value << "\t" ;
00048     }
00049   cout << endl ;
00050 } ;
00051 //___________________
00052 dProjector::~dProjector() 
00053 { 
00054   delete [] mFitParameters ;
00055 } ;
00056 //___________________
00057 TH1D* dProjector::get1dProjection(TString axis, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax)
00058 {
00059   // get ranges of the ratio 3d histogram 
00060   int xminBin = mRatio->GetXaxis()->FindFixBin(xmin) ;
00061   int xmaxBin = mRatio->GetXaxis()->FindFixBin(xmax) ;
00062   int yminBin = mRatio->GetYaxis()->FindFixBin(ymin) ;
00063   int ymaxBin = mRatio->GetYaxis()->FindFixBin(ymax) ;
00064   int zminBin = mRatio->GetZaxis()->FindFixBin(zmin) ;
00065   int zmaxBin = mRatio->GetZaxis()->FindFixBin(zmax) ;
00066 
00067   // get dimensions of the projection histos
00068   int proNbins = 0 ;
00069   double proMin = 0.0 ;
00070   double proMax = 0.0 ;
00071   char* axistitle = new char[20] ;
00072   if (axis == "x")        {  proNbins = __ABS__(xmaxBin-xminBin)+1 ;  proMin = xmin ; proMax = xmax ; strcpy(axistitle,mRatio->GetXaxis()->GetTitle()) ; }
00073   else if (axis == "y")   {  proNbins = __ABS__(ymaxBin-yminBin)+1 ;  proMin = ymin ; proMax = ymax ; strcpy(axistitle,mRatio->GetYaxis()->GetTitle()) ; }
00074   else if (axis == "z")   {  proNbins = __ABS__(zmaxBin-zminBin)+1 ;  proMin = zmin ; proMax = zmax ; strcpy(axistitle,mRatio->GetZaxis()->GetTitle()) ; }
00075   else { cout << "Wrong axis in projector, this should not happen ! "<< endl ; return 0; } ;
00076 
00077   // create 1d histo
00078   int ran = (int) (xmin*ymin*zmin*10000) ; // to distinguish between objects give them a special random id ...
00079   m1dpro = new TH1D("my1dpro","my1dpro",proNbins,proMin,proMax) ; 
00080   m1dpro->SetXTitle(axistitle) ;
00081   m1dpro->SetTitle(axistitle) ;
00082   char restitle[20] ; sprintf(restitle,"fit%s%d\n",axis.Data(),ran) ;
00083   m1dpro->SetName(restitle) ;
00084   m1dfit = new TH1D("my1dfit","my1dfit",proNbins,proMin,proMax) ;
00085   char fittitle[20] ; sprintf(fittitle,"fit%s%d\n",axis.Data(),ran) ;
00086   m1dfit->SetName(fittitle);
00087   m1dfit->SetLineColor(2) ;
00088   m1dfit->SetLineStyle(2) ;
00089   m1dfit->SetLineWidth(5) ;
00090   m1dnor = new TH1D("my1dnor","my1dnor",proNbins,proMin,proMax) ;
00091   char normtitle[20] ; sprintf(normtitle,"norm%s%d\n",axis.Data(),ran) ;
00092   m1dnor->SetName(normtitle);
00093 
00094   // fill histo
00095   for(int indexX = xminBin ; indexX <= xmaxBin ; indexX++)
00096     {
00097       for(int indexY = yminBin ; indexY <= ymaxBin ; indexY++)
00098         {
00099           for(int indexZ = zminBin ; indexZ <= zmaxBin ; indexZ++)
00100             {
00101               // get the content of the 3d cell
00102               double cellcontent = mRatio->GetBinContent(indexX,indexY,indexZ) ;
00103               // get the content of the 3d cell using the fit function
00104               TVector3 position(mRatio->GetXaxis()->GetBinCenter(indexX),
00105                                 mRatio->GetYaxis()->GetBinCenter(indexY),
00106                                 mRatio->GetZaxis()->GetBinCenter(indexZ));
00107               double cellContentFromFit = mFitter->ykpCorrelationFunction(position,mFitParameters) ;
00108               
00109               
00110               // fill the 1d histos
00111               int proBin = - 100 ;
00112               if (axis == "x")        { proBin = m1dpro->GetXaxis()->FindFixBin(mRatio->GetXaxis()->GetBinCenter(indexX)); }  
00113               else if (axis == "y")   { proBin = m1dpro->GetXaxis()->FindFixBin(mRatio->GetYaxis()->GetBinCenter(indexY)); }
00114               else if (axis == "z")   { proBin = m1dpro->GetXaxis()->FindFixBin(mRatio->GetZaxis()->GetBinCenter(indexZ)); }
00115               
00116               if(cellcontent>0.5 && cellcontent <1.5)
00117                 {
00118                   m1dpro->AddBinContent(proBin,cellcontent) ;
00119                   m1dfit->AddBinContent(proBin,cellContentFromFit) ;
00120                   
00121                   if (axis == "z" && proBin<10) cout << proBin << " " ;
00122                   // count the number of fillings to normalize the projection to 1
00123                   m1dnor->AddBinContent(proBin,1.0) ;
00124                 }
00125             }
00126         }
00127     }
00128   // normalize fit and projection histos
00129   m1dpro->Divide(m1dnor) ;
00130   m1dfit->Divide(m1dnor) ;
00131 
00132   // return projection
00133   return m1dpro ;
00134 }
00135 //___________________
00136 TH2D* dProjector::get2dProjection(TString axis, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax)
00137 {
00138     // get ranges of the ratio 3d histogram 
00139     int xminBin = mRatio->GetXaxis()->FindFixBin(xmin) ;
00140     int xmaxBin = mRatio->GetXaxis()->FindFixBin(xmax) ;
00141     int yminBin = mRatio->GetYaxis()->FindFixBin(ymin) ;
00142     int ymaxBin = mRatio->GetYaxis()->FindFixBin(ymax) ;
00143     int zminBin = mRatio->GetZaxis()->FindFixBin(zmin) ;
00144     int zmaxBin = mRatio->GetZaxis()->FindFixBin(zmax) ;
00145     
00146     // get dimensions of the projection histos
00147     int proNbinsX = 0 ;
00148     double proMinX = 0.0 ;
00149     double proMaxX = 0.0 ;
00150     char* axistitleX = new char[20] ;
00151     int proNbinsY = 0 ;
00152     double proMinY = 0.0 ;
00153     double proMaxY = 0.0 ;
00154     char* axistitleY = new char[20] ;
00155 
00156     if (axis == "x")        
00157         {  
00158             proNbinsX = __ABS__(ymaxBin-yminBin)+1 ;  proMinX = ymin ; proMaxX = ymax ; strcpy(axistitleX,mRatio->GetYaxis()->GetTitle()) ; 
00159             proNbinsY = __ABS__(zmaxBin-zminBin)+1 ;  proMinY = zmin ; proMaxY = zmax ; strcpy(axistitleY,mRatio->GetZaxis()->GetTitle()) ; 
00160         }
00161     else if (axis == "y")   
00162         {  
00163             proNbinsX = __ABS__(zmaxBin-zminBin)+1 ;  proMinX = zmin ; proMaxX = zmax ; strcpy(axistitleX,mRatio->GetZaxis()->GetTitle()) ; 
00164             proNbinsY = __ABS__(xmaxBin-xminBin)+1 ;  proMinY = xmin ; proMaxY = xmax ; strcpy(axistitleY,mRatio->GetXaxis()->GetTitle()) ;
00165         }
00166     else if (axis == "z")   
00167         {  
00168             proNbinsX = __ABS__(xmaxBin-xminBin)+1 ;  proMinX = xmin ; proMaxX = xmax ; strcpy(axistitleX,mRatio->GetXaxis()->GetTitle()) ; 
00169             proNbinsY = __ABS__(ymaxBin-yminBin)+1 ;  proMinY = ymin ; proMaxY = ymax ; strcpy(axistitleY,mRatio->GetYaxis()->GetTitle()) ;
00170         }
00171 
00172     else { cout << "Wrong axis in projector, this should not happen ! "<< endl ; return 0; } ;
00173     
00174     // create 2d histo
00175     int ran = (int) (xmin*ymin*zmin*10000) ; // to distinguish between objects give them a special random id ...
00176     // projection
00177     m2dpro = new TH2D("my2dpro","my2dpro",proNbinsX,proMinX,proMaxX,proNbinsY,proMinY,proMaxY) ; 
00178     char resName[20] ; sprintf(resName,"res%s%d\n",axis.Data(),ran) ;  m2dpro->SetName(resName) ;
00179     m2dpro->SetXTitle(axistitleX) ;
00180     m2dpro->SetYTitle(axistitleY) ;
00181     char restitle[20] ; sprintf(restitle,"%s\n",axis.Data()) ;m2dpro->SetTitle(restitle) ;
00182     // fit
00183     m2dfit = new TH2D("my2dfit","my2dfit",proNbinsX,proMinX,proMaxX,proNbinsY,proMinY,proMaxY) ;
00184     char fittitle[20] ; sprintf(fittitle,"fit%s%d\n",axis.Data(),ran) ; m2dfit->SetName(fittitle) ;
00185     // norm
00186     m2dnor = new TH2D("my2dnor","my2dnor",proNbinsX,proMinX,proMaxX,proNbinsY,proMinY,proMaxY) ;
00187     char normtitle[20] ; sprintf(normtitle,"norm%s%d\n",axis.Data(),ran) ; m2dnor->SetName(normtitle) ;
00188     
00189     // fill histo
00190     for(int indexX = xminBin ; indexX <= xmaxBin ; indexX++)
00191         {
00192             for(int indexY = yminBin ; indexY <= ymaxBin ; indexY++)
00193                 {
00194                     for(int indexZ = zminBin ; indexZ <= zmaxBin ; indexZ++)
00195                         {
00196                             // get the content of the 3d cell
00197                             double cellcontent = mRatio->GetBinContent(indexX,indexY,indexZ) ;
00198                             // get the content of the 3d cell using the fit function
00199                             TVector3 position(mRatio->GetXaxis()->GetBinCenter(indexX),
00200                                               mRatio->GetYaxis()->GetBinCenter(indexY),
00201                                               mRatio->GetZaxis()->GetBinCenter(indexZ));
00202                             double cellContentFromFit = mFitter->ykpCorrelationFunction(position,mFitParameters) ;
00203                             
00204                             
00205                             // fill the 1d histos
00206                             double proX = 100000.0 ;
00207                             double proY = 100000.0 ;
00208                             if (axis == "x")        
00209                                 { 
00210                                     proX = mRatio->GetYaxis()->GetBinCenter(indexY) ; 
00211                                     proY = mRatio->GetZaxis()->GetBinCenter(indexZ) ;
00212                                 }  
00213                             else if (axis == "y")   
00214                                 { 
00215                                     proX = mRatio->GetZaxis()->GetBinCenter(indexZ) ; 
00216                                     proY = mRatio->GetXaxis()->GetBinCenter(indexX) ;
00217                                 }
00218                             else if (axis == "z")   
00219                                 { 
00220                                     proX = mRatio->GetXaxis()->GetBinCenter(indexX) ; 
00221                                     proY = mRatio->GetYaxis()->GetBinCenter(indexY) ;
00222                                 }
00223                             
00224                             if(cellcontent>0.5 && cellcontent <1.5)
00225                                 {
00226                                     // now fill 2d histos 
00227                                     m2dpro->Fill(proX,proY,cellcontent) ;
00228                                     m2dfit->Fill(proX,proY,cellContentFromFit) ;
00229                                     // count the number of fillings to normalize the projection to 1
00230                                     m2dnor->Fill(proX,proY,1.0) ;
00231                                 }
00232                         }
00233                 }
00234         }
00235     // normalize fit and projection histos
00236     m2dpro->Divide(m2dnor) ;
00237     m2dfit->Divide(m2dnor) ;
00238     
00239     // return projection
00240     return m2dpro ;
00241     
00242 }