StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFcsPointMaker.cxx
1 // $Id: StFcsPointMaker.cxx,v 1.1 2021/03/30 13:40:10 akio Exp $
2 
3 #include "StFcsPointMaker.h"
4 #include "StLorentzVectorF.hh"
5 
6 #include "StMessMgr.h"
7 #include "StEvent/StEventTypes.h"
8 #include "StEvent/StFcsHit.h"
9 #include "StEvent/StFcsCluster.h"
10 #include "StEvent/StFcsPoint.h"
11 #include "StFcsDbMaker/StFcsDb.h"
12 
13 #include "StMuDSTMaker/COMMON/StMuTypes.hh"
14 #include "StMuDSTMaker/COMMON/StMuDst.h"
15 #include "StMuDSTMaker/COMMON/StMuEvent.h"
16 #include "tables/St_vertexSeed_Table.h"
17 
18 #include <cmath>
19 #include <limits>
20 #include "TMath.h"
21 #include "TVector2.h"
22 
23 namespace{
24  // shower shape parameters
25  static const int NPMAX=60;
26  std::array<double,NPMAX> mShowerShapeParameters;
27 
28  // tower data for current working cluster
29  static int mNHit; // # of hit in current working cluster
30  static float mEtot; // total energy of the cluster
31  std::vector<float> mX; // x of hit
32  std::vector<float> mY; // y of hit
33  std::vector<float> mE; // E of hit
34 }
35 
36 ClassImp(StFcsPointMaker)
37 
38 StFcsPointMaker::StFcsPointMaker(const char* name) : StMaker(name) , mMinuit(7) {
39  mMinuit.SetPrintLevel(-1);
40 }
41 
42 StFcsPointMaker::~StFcsPointMaker() { }
43 
44 void StFcsPointMaker::Clear(Option_t* option) {
45  StMaker::Clear(option);
46 }
47 
48 int StFcsPointMaker::InitRun(int runNumber) {
49  // Ensure we can access database information
50  LOG_DEBUG << "StFcsPointMaker initializing run" << endm;
51  mDb = static_cast<StFcsDb*>(GetDataSet("fcsDb"));
52  if (!mDb) {
53  LOG_ERROR << "StFcsPointMaker initializing failed due to no StFcsDb" << endm;
54  return kStErr;
55  }
56  return StMaker::InitRun(runNumber);
57 }
58 
59 void StFcsPointMaker::setShowerShapeParameters(int det){
60  double unused=0.0;
61  double width = mDb->getXWidth(det);
62  double scl=mShowerShapeScale; //scale shower shape to cell width and adjuistable scale factor
63  if(mShowerShape==0){
64  //Original single slice shower shap
65  scl *= 1.0/3.849; //optimum from FMS study, scale it to cell size
66  mShowerShapeParameters= { width, 1.0708, 0.167773, -0.238578, 0.535845*scl, 0.850233*scl, 2.38264*scl, 1.0, unused, unused,
67  unused, unused, unused, unused, unused, unused, unused, unused, unused, unused,
68  unused, unused, unused, unused, unused, unused, unused, unused, unused, unused,
69  unused, unused, unused, unused, unused, unused, unused, unused, unused, unused,
70  unused, unused, unused, unused, unused, unused, unused, unused, unused, unused,
71  unused, unused, unused, unused, unused, unused, unused, unused, unused, unused};
72  }else if(mShowerShape==1){
73  //Yuxi's 6 slices * 3 gaus
74  scl *= 0.8/3.849; //optimum from FMS study, scale it to cell size
75  double a11=0.998438; double a12=0.222782; double a13=-0.22122; double b11=0.177028;double b12=0.000473222;double b13=0.178897;double w1 = 0.0372556;
76  double a21=1.07711 ; double a22=-0.0281385; double a23=-0.0489747; double b21=0.199964;double b22=3.5021; double b23=2.35246; double w2 = 0.202699;
77  double a31=1.07901 ; double a32=0.0650143; double a33=-0.144025; double b31=0.446845;double b32=0.00544512; double b33=1.64565; double w3 = 0.293878;
78  double a41=0.922174; double a42=0.0778254; double a43=1.07474e-07;double b41=0.593804;double b42=0.6199; double b43=3.49798; double w4 = 0.236854;
79  double a51=0.999849; double a52=0.000151185;double a53=2.20244e-07;double b51=0.949953;double b52=1.84451; double b53=3.40799; double w5 = 0.146041;
80  double a61=0.997454; double a62=0.00254497; double a63=1.02127e-06;double b61=1.43387; double b62=2.91155; double b63=3.4484; double w6 = 0.0832717;
81  mShowerShapeParameters= {width, a11*w1, a12*w1, a13*w1, b11*scl, b12*scl, b13*scl, unused ,unused, unused,
82  width, a21*w2, a22*w2, a23*w2, b21*scl, b22*scl, b23*scl, unused ,unused, unused,
83  width, a31*w3, a32*w3, a33*w3, b31*scl, b32*scl, b33*scl, unused ,unused, unused,
84  width, a41*w4, a42*w4, a43*w4, b41*scl, b42*scl, b43*scl, unused ,unused, unused,
85  width, a51*w5, a52*w5, a53*w5, b51*scl, b52*scl, b53*scl, unused ,unused, unused,
86  width, a61*w6, a62*w6, a63*w6, b61*scl, b62*scl, b63*scl, unused ,unused, unused};
87  }else if(mShowerShape==2){
88  //Zhanwen's small 45GeV new (6 slices * 2 gaus)
89  scl *= 0.6/3.849; //optimum from FMS study, scale it to cell size
90  double a1S[6]={0.0303644,0.212191,0.277429,0.0370035,0.0524404,0.00844062};
91  double a2S[6]={0.00122867,0.105355,0.10538,0.152656,0.00664331,0.0108688};
92  double b1S[6]={0.403493,0.514546,0.672826,1.82344,0.727991,1.48785};
93  double b2S[6]={0.270492,0.514593,0.672655,0.644871,4.32003,0.25};
94  mShowerShapeParameters= {width, a1S[0], a2S[0], 0.0, b1S[0]*scl, b2S[0]*scl, 0.0, unused, unused, unused,
95  width, a1S[1], a2S[1], 0.0, b1S[1]*scl, b2S[1]*scl, 0.0, unused, unused, unused,
96  width, a1S[2], a2S[2], 0.0, b1S[2]*scl, b2S[2]*scl, 0.0, unused, unused, unused,
97  width, a1S[3], a2S[3], 0.0, b1S[3]*scl, b2S[3]*scl, 0.0, unused, unused, unused,
98  width, a1S[4], a2S[4], 0.0, b1S[4]*scl, b2S[4]*scl, 0.0, unused, unused, unused,
99  width, a1S[5], a2S[5], 0.0, b1S[5]*scl, b2S[5]*scl, 0.0, unused, unused, unused};
100  }else if(mShowerShape==3){
101  //Zhanwen's large 45GeV new (6 slices & 2 gaus)
102  scl *= 0.8/5.812; //optimum from FMS study, scale it to cell size
103  double a1L[6]={0.0275364,0.200363,0.277157,0.0391611,0.0590757,0.0101089};
104  double a2L[6]={0.000429808,0.0991777,0.104781,0.161916,0.00764026,0.012653};
105  double b1L[6]={0.515974,0.661722,0.865167,2.35237,0.932038,1.87933};
106  double b2L[6]={0.53531,0.661519,0.865226,0.828017,5.49041,0.321139};
107  mShowerShapeParameters= {width, a1L[0], a2L[0], 0.0, b1L[0]*scl, b2L[0]*scl, 0.0, unused, unused, unused,
108  width, a1L[1], a2L[1], 0.0, b1L[1]*scl, b2L[1]*scl, 0.0, unused, unused, unused,
109  width, a1L[2], a2L[2], 0.0, b1L[2]*scl, b2L[2]*scl, 0.0, unused, unused, unused,
110  width, a1L[3], a2L[3], 0.0, b1L[3]*scl, b2L[3]*scl, 0.0, unused, unused, unused,
111  width, a1L[4], a2L[4], 0.0, b1L[4]*scl, b2L[4]*scl, 0.0, unused, unused, unused,
112  width, a1L[5], a2L[5], 0.0, b1L[5]*scl, b2L[5]*scl, 0.0, unused, unused, unused};
113  }
114  if(mShowerShape>0){ //get slope factor from vertex=(0,0,0) to z of 6 slices
115  StThreeVectorD off = mDb->getDetectorOffset(det);
116  double dz=mDb->getZDepth(det);
117  double smax=mDb->getShowerMaxZ(det);
118  double a = mDb->getDetectorAngle(det) / 180.0 * M_PI;
119  double z0 = off.z()/cos(a);
120  for(int i=0; i<6; i++){
121  double z = (i+0.5)*dz/6.0;
122  mShowerShapeParameters[i*10+7]=(z0 + z)/(z0+smax);
123  }
124  }
125  if(GetDebug()>0) {
126  for(int i=0; i<6; i++){
127  LOG_INFO << Form("Shower Shape Parameters det=%1d slice=%1d : ",det,i);
128  for(int j=0; j<10; j++) {LOG_INFO << Form(" %10.6f",mShowerShapeParameters[i*10+j]);}
129  LOG_INFO << endm;
130  }
131  }
132 }
133 
134 int StFcsPointMaker::Make() {
135  LOG_DEBUG << "StFcsPointMaker Make!!!" << endm;
136 
137  StEvent* event = static_cast<StEvent*>(GetInputDS("StEvent"));
138  mFcsCollection=0;
139  if (event) mFcsCollection = event->fcsCollection();
140  if(!mFcsCollection) {
141  LOG_WARN << "StFcsPointMaker did not find fcsCollection in StEvent" << endm;
142  return kStWarn;
143  }
144 
145  for(int det=0; det<=kFcsEcalSouthDetId; det++) {
146  fitClusters(det);
147  }
148  if(GetDebug()>0) mFcsCollection->print(3);
149  return kStOk;
150 }
151 
152 void StFcsPointMaker::fitClusters(int det) {
153  StSPtrVecFcsCluster& clusters = mFcsCollection->clusters(det);
154  StSPtrVecFcsPoint& points = mFcsCollection->points(det);
155 
156  points.clear(); //clear all points
157  int nclu=clusters.size();
158  for(int i=0; i<nclu; i++) clusters[i]->points().clear(); //reset point pointer from cluster
159 
160  setShowerShapeParameters(det);
161 
162  StFcsPoint point0,point1,point2;
163  for(int i=0; i<nclu; i++){ //loop over all clusters
164  StFcsCluster* c=clusters[i];
165  mNHit=c->nTowers();
166  mEtot=c->energy();
167  //mX=new float[mNHit];
168  //mY=new float[mNHit];
169  //mE=new float[mNHit];
170  mE.clear(); mX.clear(); mY.clear();
171  for(int j=0; j<mNHit; j++){
172  StFcsHit* h=c->hits()[j];
173  float x,y;
174  mDb->getLocalXYinCell(h, x, y);
175  mE.push_back(h->energy());
176  mX.push_back(x);
177  mY.push_back(y);
178  //mE[j]=h->energy();
179  //mX[j]=x;
180  //mY[j]=y;
181  }
182  double chi1=std::numeric_limits<double>::max();
183  double chi2=std::numeric_limits<double>::max();
184  switch(c->category()){
185  case 0:
186  chi1 = fit1PhotonCluster(c,&point0);
187  chi2 = fit2PhotonCluster(c,&point1,&point2);
188  break;
189  case 1:
190  chi1 = fit1PhotonCluster(c,&point0);
191  break;
192  case 2:
193  chi2 = fit2PhotonCluster(c,&point1,&point2);
194  break;
195  default:
196  chi1 = chi2 = 0;
197  LOG_WARN << "Unknown cluster category=" << c->category() << endm;
198  }
199  // sotre chi2 for both
200  c->setChi2Ndf1Photon(chi1);
201  c->setChi2Ndf2Photon(chi2);
202  // pick better one
203  if(chi1<chi2){
204  StFcsPoint* p0=new StFcsPoint(point0);
205  p0->setDetectorId(det);
206  p0->setCluster(c);
207  p0->setNParentClusterPhotons(1);
208  StThreeVectorD xyz=mDb->getStarXYZfromColumnRow(det,p0->x(),p0->y());
209  p0->setXYZ(xyz);
210  p0->setFourMomentum(mDb->getLorentzVector(xyz,p0->energy()));
211  mFcsCollection->addPoint(det,p0);
212  c->addPoint(p0);
213  }else{
214  StFcsPoint* p1=new StFcsPoint(point1);
215  p1->setDetectorId(det);
216  p1->setCluster(c);
217  p1->setNParentClusterPhotons(2);
218  StThreeVectorD xyz1=mDb->getStarXYZfromColumnRow(det,p1->x(),p1->y());
219  p1->setXYZ(xyz1);
220  p1->setFourMomentum(mDb->getLorentzVector(xyz1,p1->energy()));
221  mFcsCollection->addPoint(det,p1);
222 
223  StFcsPoint* p2=new StFcsPoint(point2);
224  p2->setDetectorId(det);
225  p2->setCluster(c);
226  p2->setNParentClusterPhotons(2);
227  StThreeVectorD xyz2=mDb->getStarXYZfromColumnRow(det,p2->x(),p2->y());
228  p2->setXYZ(xyz2);
229  p2->setFourMomentum(mDb->getLorentzVector(xyz2,p2->energy()));
230  mFcsCollection->addPoint(det,p2);
231 
232  c->addPoint(p1,p2);
233  }
234  }
235 
236  //loop over all found points and fill fourMomentum
237  int np = points.size();
238  for(int i=0; i<np; i++){
239  StFcsPoint* p=points[i];
240  StThreeVectorF xyz = mDb->getStarXYZfromColumnRow(det,p->x(),p->y());
241  p->setFourMomentum(mDb->getLorentzVector(xyz,p->energy(),0.0));
242  }
243 }
244 
245 
246 double StFcsPointMaker::fit1PhotonCluster(StFcsCluster* c, StFcsPoint* p){
247  int err;
248  double chi2 = -1.0;
249  mMinuit.SetFCN(minimizationFunctionNPhoton);
250  mMinuit.mncler();
251  if(mMinuit.GetNumFixedPars() > 0) {mMinuit.mnfree(0);}
252  double x=c->x();
253  double y=c->y();
254  double e=c->energy();
255  double dx=m_PH1_Delta_X; //0.5
256  double de=m_PH1_Delta_E; //1.15
257  const std::vector<TString> names = {"np", "x", "y", "e"};
258  double one=1.0;
259  double zero=0.0;
260  mMinuit.DefineParameter(0, "np", one, zero, one, one);
261  mMinuit.DefineParameter(1, "x", x, dx/5, x-dx, x+dx);
262  mMinuit.DefineParameter(2, "y", y, dx/5, y-dx, y+dx);
263  mMinuit.DefineParameter(3, "e", e, e*(1-de)/5, e/de, e*de);
264  //double fval; mMinuit.mnprin(1,fval);
265  //mMinuit.FixParameter(0);
266  if(m_PH1_FixEnergy==1) mMinuit.FixParameter(3);
267 
268  //run minimization
269  err = -1;
270  double arguments[2]={1000.0,1.0}; //Max calls and tolerance
271  mMinuit.mnexcm("MIGRAD", arguments, 2, err);
272 
273  if(mMinuit.GetStatus() == 0){
274  double par[4],perr[4],g[4];
275  for(int i=0; i<4; i++){ mMinuit.GetParameter(i, par[i], perr[i]); }
276  p->setX(par[1]);
277  p->setY(par[2]);
278  p->setEnergy(par[3]);
279  mMinuit.Eval(4, g, chi2, par, err);
280  }
281  return chi2;
282 }
283 
284 double StFcsPointMaker::fit2PhotonCluster(StFcsCluster* c, StFcsPoint* p1, StFcsPoint* p2){
285  int err;
286  double chi2 = -1.0;
287  mMinuit.SetFCN(minimizationFunction2Photon);
288  mMinuit.mncler();
289  if(mMinuit.GetNumFixedPars() > 0) {mMinuit.mnfree(0);}
290  double x=c->x();
291  double y=c->y();
292  double smax= c->sigmaMax();
293  // double smin= c->sigmaMin();
294  double d=std::max(m_PH2_StartDggFactor*smax,0.5); //1.1 = 2.2/2.0
295  double t=c->theta();
296  double z=0.0;
297  double e=c->energy();
298  double dx=m_PH2_Delta_X; //0.2
299  double ddl=d*m_PH2_Low_Dgg; //0.8
300  double ddh=d*m_PH2_High_Dgg; //3.0
301  double dt=m_PH2_MaxTheta_F; //TMath::PiOver2()
302  double de=m_PH2_Delta_E; //1.05
303  mMinuit.mnparm(0, "np", 2.0, 0.0, 2.0, 2.0, err);
304  mMinuit.mnparm(1, "x", x, dx/5, x-dx, x+dx, err);
305  mMinuit.mnparm(2, "y", y, dx/5, y-dx, y+dx, err);
306  mMinuit.mnparm(3, "dgg", d, d*0.05, ddl, ddh, err);
307  mMinuit.mnparm(4, "theta", t, dt/5, t-dt, t+dt, err);
308  mMinuit.mnparm(5, "zgg", z, 0.1, -0.99, 0.99, err);
309  mMinuit.mnparm(6, "etot", e, e*(1-de)/5, e/de, e*de, err);
310  mMinuit.FixParameter(0);
311  // Fix E_total and theta, we don't want these to be free parameters
312  if(m_PH2_FixTheta==1) mMinuit.FixParameter(4);
313  if(m_PH2_FixEnergy==1) mMinuit.FixParameter(6);
314 
315  //run minimization
316  err = -1;
317  double arguments[2]={1000.0,1.0}; //Max calls and tolerance
318  mMinuit.mnexcm("MIGRAD", arguments, 2, err);
319 
320  if(mMinuit.GetStatus() == 0) {
321  double par[7],perr[7],g[7];
322  for(int i=0; i<7; i++){ mMinuit.GetParameter(i, par[i], perr[i]); }
323  double x=par[1];
324  double y=par[2];
325  double d=par[3];
326  double t=par[4];
327  double z=par[5];
328  double e=par[6];
329  double x1 = x + cos(t)*d*(1.0-z)/2.0;
330  double y1 = y + sin(t)*d*(1.0-z)/2.0;
331  double e1 = e*(1.0+z)/2.0;
332  double x2 = x - cos(t)*d*(1.0+z)/2.0;
333  double y2 = y - sin(t)*d*(1.0+z)/2.0;
334  double e2 = e*(1.0-z)/2.0;
335  p1->setX(x1);
336  p1->setY(y1);
337  p1->setEnergy(e1);
338  p2->setX(x2);
339  p2->setY(y2);
340  p2->setEnergy(e2);
341  mMinuit.Eval(7, g, chi2, par, 1);
342  }
343  return chi2;
344 }
345 
346 // Minimization function to be called from TMinuit
347 void StFcsPointMaker::minimizationFunctionNPhoton(int& npara, double* grad, double& fval, double* para, int){
348  fval = 0.0;
349  const int nPhotons = static_cast<int>(para[0]);
350  for(int i=0; i<mNHit; i++){
351  double expected = 0.0;
352  for (int j=0; j<nPhotons; j++){
353  int k = 3*j;
354  expected += para[k+3] * energyDepositionInTower(mX[i], mY[i], para[k+1], para[k+2]);
355  }
356  const double measured = mE[i];
357  const double deviation = measured - expected;
358  const double ratio = measured / mEtot;
359  const double err = 0.01 + 0.03 * pow(ratio, 1.0-0.001*mEtot) * pow(1.0-ratio, 1.0-0.007*mEtot)*mEtot;
360  fval += deviation * deviation / err;
361  //LOG_INFO << Form("i=%2d M=%6.2f E=%6.2f D=%6.2f R=%6.2f E=%6.2f F=%6.2f",i,measured,expected,deviation,ratio,err,deviation * deviation / err)<<endm;
362  }
363  fval = std::max(fval, 0.);
364  //LOG_INFO << Form("FVAL=%6.2f",fval)<<endm;
365 }
366 
367 // Minimization function to be called from TMinuit for 2 photon case only
368 void StFcsPointMaker::minimizationFunction2Photon(int& npara, double* grad, double& fval, double* para, int){
369  // Only need to translate into the old parameterization
370  const double dgg = para[3];
371  const double zgg = para[5];
372  const double angle = para[4];
373  double oldPara[7] = {
374  para[0], // Number of photons, unchanged
375  para[1] + cos(angle) * dgg * (1 - zgg) / 2.0, // x 1
376  para[2] + sin(angle) * dgg * (1 - zgg) / 2.0, // y 1
377  para[6] * (1 + zgg) / 2.0, // e 1
378  para[1] - cos(angle) * dgg * (1 + zgg) / 2.0, // x 2
379  para[2] - sin(angle) * dgg * (1 + zgg) / 2.0, // y 2
380  para[6] * (1 - zgg) / 2.0 // e 2
381  };
382  // Now call the regular minimization function with the translated parameters
383  minimizationFunctionNPhoton(npara, grad, fval, oldPara, 0);
384 }
385 
386 // getting energy deposition in a tower by summing up 6 z slices
387 // x,y are cell center position
388 // xun,yun are photon position in cell coordinate at Z=shower max
389 // xc,yc are photon position at z of the slice
390 // (scale factor = mShowerShapeParameters[istart+7] calculated in setShowerShapeParameters()
391 double StFcsPointMaker::energyDepositionInTower(double x, double y,double xun, double yun){
392  double sum = 0.0;
393  for(int i=0; i<6; i++){
394  int istart = i*10;
395  if(mShowerShapeParameters[istart]>0.0){
396  double xc = xun * mShowerShapeParameters[istart+7];
397  double yc = yun * mShowerShapeParameters[istart+7];
398  sum += energyDepositionInTowerSingleLayer(x-xc,y-yc,&mShowerShapeParameters[istart]);
399  }
400  }
401  return sum;
402 }
403 
StFcsDb::getLorentzVector
StLorentzVectorD getLorentzVector(const StThreeVectorD &xyz, float energy, float zVertex=0.0)
Get get 4 vector assuing m=0 and taking beamline from DB.
Definition: StFcsDb.cxx:895
StMaker::Clear
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StMaker.cxx:634
StFcsDb::getDetectorOffset
StThreeVectorD getDetectorOffset(int det, double zdepth=-1) const
Utility functions related to DetectorPosition.
Definition: StFcsDb.cxx:578
StFcsDb::getStarXYZfromColumnRow
StThreeVectorD getStarXYZfromColumnRow(int det, float col, float row, float FcsZ=-1.0) const
get the STAR frame cooridnates from other way
Definition: StFcsDb.cxx:863
StFcsDb::getZDepth
float getZDepth(int det) const
get the Y width of the cell
Definition: StFcsDb.cxx:639
StFcsDb::getXWidth
float getXWidth(int det) const
get the angle of the detector
Definition: StFcsDb.cxx:627
StFcsDb::getLocalXYinCell
void getLocalXYinCell(StFcsHit *hit, float &x, float &y) const
getting XY in local cell coordinate
Definition: StFcsDb.cxx:652
StFcsPointMaker::Clear
void Clear(Option_t *option="")
User defined functions.
Definition: StFcsPointMaker.cxx:44
StFcsDb::getDetectorAngle
float getDetectorAngle(int det) const
This is the vector normal to the detector plane.
Definition: StFcsDb.cxx:617
kStWarn
Definition: Stypes.h:42
kStErr
Definition: Stypes.h:44
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5