StRoot: StGammaMaker/StGammaFitter.cxx Source File

00001 //
00002 // Pibero Djawotho <pibero@indiana.edu>
00003 // Indiana University
00004 // 28 July 2007
00005 //
00006 
00007 // ROOT
00008 #include "TFile.h"
00009 #include "TH1.h"
00010 #include "TF1.h"
00011 #include "TCanvas.h"
00012 
00013 // STAR
00014 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00015 #include "StEEmcUtil/EEmcSmdMap/EEmcSmdMap.h"
00016 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
00017 #include "StEmcUtil/geometry/StEmcGeom.h"
00018 
00019 // Local
00020 #include "StGammaEvent.h"
00021 #include "StGammaCandidate.h"
00022 #include "StGammaFitterResult.h"
00023 #include "StGammaFitter.h"
00024 
00025 ClassImp(StGammaFitter);
00026 
00027 StGammaFitter* StGammaFitter::mInstance = 0;
00028 TH1* StGammaFitter::hU = 0;
00029 TH1* StGammaFitter::hV = 0;
00030 TF1* StGammaFitter::fFit[2];
00031 TF1* StGammaFitter::fResidualCut = 0;
00032 int StGammaFitter::mNdf = 0;
00033 TCanvas* StGammaFitter::mCanvas = 0;
00034 TF1* StGammaFitter::mShowerShapes[3];
00035 
00036 StGammaFitter::StGammaFitter()
00037 {
00038   hU = new TH1F("hU", "hU", 288, 0, 288);
00039   hV = new TH1F("hV", "hV", 288, 0, 288);
00040 
00041   fResidualCut = new TF1("fResidualCut", "pol2", 0, 200);
00042   fResidualCut->SetParameters(100, 0, 0.05);
00043 
00044   // Initialize shower shapes with fits from Ilya based on his gamma-jet analysis for
00045   // different preshower conditions.
00046   // See http://www.star.bnl.gov/protected/spin/seluzhen/gammaJet/2008.05.27/showerShapes_comparison.html
00047   static const char* formula[] = {
00048     "[0]*(0.669864*exp(-0.5*sq((x-[1])/0.574864))+0.272997*exp(-0.5*sq((x-[1])/-1.84608))+0.0585682*exp(-0.5*sq((x-[1])/5.49802)))", // pre1=0, pre2=0
00049     "[0]*(0.0694729*exp(-0.5*sq((x-[1])/5.65413))+0.615724*exp(-0.5*sq((x-[1])/0.590723))+0.314777*exp(-0.5*sq((x-[1])/2.00192)))", // pre1=0, pre2>0
00050     "[0]*(0.0955638*exp(-0.5*sq((x-[1])/5.59675))+0.558661*exp(-0.5*sq((x-[1])/0.567596))+0.345896*exp(-0.5*sq((x-[1])/1.9914)))" // pre1>0, pre2>0
00051   };
00052 
00053   for (int i = 0; i < 3; ++i) {
00054     const char* name = Form("mShowerShapes%d", i);
00055     mShowerShapes[i] = new TF1(name, formula[i], 0, 288);
00056   }
00057 }
00058 
00059 StGammaFitter::~StGammaFitter()
00060 {
00061   // Never called???
00062   delete hU;
00063   delete hV;
00064   delete fResidualCut;
00065 }
00066 
00067 StGammaFitter* StGammaFitter::instance()
00068 {
00069   if (!mInstance) mInstance = new StGammaFitter;
00070   return mInstance;
00071 }
00072 
00073 int StGammaFitter::fit(StGammaCandidate* candidate, StGammaFitterResult* fits, Int_t plane)
00074 {
00075   // Require at least 5 strips in each SMD plane
00076   if (candidate->numberOfSmdu() < 5) return 9;
00077   if (candidate->numberOfSmdv() < 5) return 9;
00078 
00079   // Clear fit results
00080   //memset(fits, 0, 1 * sizeof(StGammaFitterResult));
00081 
00082   TF1* fit = 0;
00083 
00084   if (candidate->pre1Energy() == 0 && candidate->pre2Energy() == 0)
00085     fit = mShowerShapes[0];
00086   else if (candidate->pre1Energy() == 0 && candidate->pre2Energy() > 0)
00087     fit = mShowerShapes[1];
00088   else if (candidate->pre1Energy() > 0 && candidate->pre2Energy() > 0)
00089     fit = mShowerShapes[2];
00090   else
00091     return 9;
00092 
00093   // Loop over planes
00094   //  for (int plane = 0; plane < 2; ++plane) {
00095   {
00096     static TH1* hStrips = new TH1F("hStrips", "hStrips", 288, 0, 288);
00097     hStrips->Reset();
00098 
00099     switch (plane) {
00100     case 0:
00101       for (int i = 0; i < candidate->numberOfSmdu(); ++i) {
00102         StGammaStrip* strip = candidate->smdu(i);
00103         hStrips->Fill(strip->index, strip->energy);
00104       }
00105       break;
00106     case 1:
00107       for (int i = 0; i < candidate->numberOfSmdv(); ++i) {
00108         StGammaStrip* strip = candidate->smdv(i);
00109         hStrips->Fill(strip->index, strip->energy);
00110       }
00111       break;
00112     }
00113 
00114     // Find maximum strip
00115     int mean = hStrips->GetMaximumBin();
00116 
00117     // Integrate yield from +/- 2 strips of max strip
00118     int bin1 = max(1,   mean - 2);
00119     int bin2 = min(288, mean + 2);
00120 
00121     // Fit to Ilya's shower shape
00122     float yield = hStrips->Integral(bin1, bin2);
00123     fit->SetParameters(yield, hStrips->GetBinLowEdge(mean));
00124     hStrips->Fit(fit, "WWQ", "", hStrips->GetBinLowEdge(bin1), hStrips->GetBinLowEdge(bin2));
00125 
00126     // Save fit results
00127     fits[plane].yield = yield;
00128     fits[plane].yieldError = 0;
00129     fits[plane].centroid = fit->GetParameter(1);
00130     fits[plane].centroidError = fit->GetParError(1);
00131     fits[plane].residual = residual(hStrips, fit);
00132     fits[plane].mean = hStrips->GetMean();
00133     fits[plane].rms = hStrips->GetRMS();
00134     fits[plane].chiSquare = fit->GetChisquare();
00135     fits[plane].ndf = fit->GetNDF();
00136     fits[plane].prob = fit->GetProb();
00137   }
00138 
00139   return 0;
00140 }
00141 
00142 void StGammaFitter::estimateYieldMean(TH1* h1, float& yield, float& mean)
00143 {
00144   int bin = h1->GetMaximumBin();
00145   int xmin = max(bin - 2, 1);
00146   int xmax = min(bin + 2, 288);
00147   yield = h1->Integral(xmin, xmax);
00148   mean = h1->GetBinCenter(bin);
00149 }
00150 
00151 float StGammaFitter::residual(TH1* h1, TF1* f1)
00152 {
00153   int mean = h1->FindBin(f1->GetParameter(1));
00154 
00155   int leftMin = 1;
00156   int leftMax = max(mean - 3, 1);
00157   float leftResidual = 0;
00158 
00159   for (int bin = leftMin; bin <= leftMax; ++bin) {
00160     float data = h1->GetBinContent(bin);
00161     float x = h1->GetBinCenter(bin);
00162     float fit = f1->Eval(x);
00163     leftResidual += data - fit;
00164   }
00165 
00166   int rightMin = min(mean + 3, 288);
00167   int rightMax = 288;
00168   float rightResidual = 0;
00169 
00170   for (int bin = rightMin; bin <= rightMax; ++bin) {
00171     float data = h1->GetBinContent(bin);
00172     float x = h1->GetBinCenter(bin);
00173     float fit = f1->Eval(x);
00174     rightResidual += data - fit;
00175   }
00176 
00177   return max(leftResidual, rightResidual);
00178 }
00179 
00180 double StGammaFitter::distanceToQuadraticCut(double x, double y)
00181 {
00182   double a = fResidualCut->GetParameter(0);
00183   double b = fResidualCut->GetParameter(2);
00184   double coef[] = { -x, 2*b*(a-y)+1, 0, 2*b*b };
00185   double r1, r2, r3;
00186   double x1 = 0;
00187 
00188   if (TMath::RootsCubic(coef, r1, r2, r3)) {
00189     // 1 real root r1 and 2 complex conjugates r2+i*r3 and r2-i*r3
00190     // Pick real root
00191     x1 = r1;
00192   }
00193   else {
00194     // 3 real roots r1, r2, r3
00195     // Pick positive root
00196     if (r1 >= 0)
00197       x1 = r1;
00198     else if (r2 >= 0)
00199       x1 = r2;
00200     else if (r3 >= 0)
00201       x1 = r3;
00202     else {
00203       LOG_ERROR << "Can't find positive root" << endm;
00204     }
00205   }
00206 
00207   double y1 = fResidualCut->Eval(x1);
00208   double d = hypot(x - x1, y - y1);
00209 
00210   return (y > fResidualCut->Eval(x)) ? d : -d;
00211 }
00212 
00213 float StGammaFitter::GetMaximum(TH1* h1, float xmin, float xmax)
00214 {
00215   int bin1 = h1->FindBin(xmin);
00216   int bin2 = h1->FindBin(xmax);
00217 
00218   float ymax = 0;
00219 
00220   for (int bin = bin1; bin <= bin2; ++bin) {
00221     float y = h1->GetBinContent(bin);
00222     if (y > ymax) ymax = y;
00223   }
00224 
00225   return ymax;
00226 }