StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFcsTrackMatchMaker.cxx
1 // \class StFcsTrackMatchMaker
2 // \author Akio Ogawa
3 //
4 // $Id: StFcsTrackMatchMaker.cxx,v 1.1 2021/03/30 13:34:15 akio Exp $
5 // $Log: StFcsTrackMatchMaker.cxx,v $
6 
7 #include "StFcsTrackMatchMaker.h"
8 #include "StEnumerations.h"
9 #include "StMessMgr.h"
10 #include "StEvent/StEvent.h"
11 #include "StEvent/StFcsCollection.h"
12 #include "StEvent/StFcsCluster.h"
13 #include "StEvent/StFwdTrackCollection.h"
14 #include "StEvent/StFwdTrack.h"
15 
16 #include "StThreeVectorF.hh"
17 #include "StFcsDbMaker/StFcsDb.h"
18 #include "StRoot/StEpdUtil/StEpdGeom.h"
19 
20 #include "TFile.h"
21 #include "TH1F.h"
22 #include "TH2F.h"
23 
24 ClassImp(StFcsTrackMatchMaker)
25 
26 StFcsTrackMatchMaker::StFcsTrackMatchMaker(const char *name) : StMaker(name)
27 {
28  mMinEnergy[0] = 0.1; // ~1/3 MIP
29  mMinEnergy[1] = 0.5; // ~1/3 MIP
30  mMaxDistance[0] = 6.0;
31  mMaxDistance[1] = 10.0;
32 }
33 
34 StFcsTrackMatchMaker::~StFcsTrackMatchMaker() {
35 
36  if (mHdx[0]) {delete mHdx[0]; mHdx[0]=nullptr;}
37  if (mHdy[0]) {delete mHdy[0]; mHdy[0]=nullptr;}
38  if (mHdr[0]) {delete mHdr[0]; mHdr[0]=nullptr;}
39  if (mNtrk[0]) {delete mNtrk[0]; mNtrk[0]=nullptr;}
40  if (mNclu[0]) {delete mNclu[0]; mNclu[0]=nullptr;}
41  if (mHdx[1]) {delete mHdx[1]; mHdx[1]=nullptr;}
42  if (mHdy[1]) {delete mHdy[1]; mHdy[1]=nullptr;}
43  if (mHdr[1]) {delete mHdr[1]; mHdr[1]=nullptr;}
44  if (mNtrk[1]) {delete mNtrk[1]; mNtrk[1]=nullptr;}
45  if (mNclu[1]) {delete mNclu[1]; mNclu[1]=nullptr;}
46  if (mNtrk[2]) {delete mNtrk[2]; mNtrk[2]=nullptr;}
47  if (mNtrk[3]) {delete mNtrk[3]; mNtrk[3]=nullptr;}
48  if (mNclu[2]) {delete mNclu[2]; mNclu[2]=nullptr;}
49  if (mNclu[3]) {delete mNclu[3]; mNclu[3]=nullptr;}
50  if (mPtEt[0]) {delete mPtEt[0]; mPtEt[0]=nullptr;}
51  if (mPtEt[1]) {delete mPtEt[1]; mPtEt[1]=nullptr;}
52  if (mPtEt2[0]) {delete mPtEt2[0]; mPtEt2[0]=nullptr;}
53  if (mPtEt2[1]) {delete mPtEt2[1]; mPtEt2[1]=nullptr;}
54  if (mCharge[0]) {delete mCharge[0]; mCharge[0]=nullptr;}
55  if (mCharge[1]) {delete mCharge[1]; mCharge[1]=nullptr;}
56  if (mCharge[2]) {delete mCharge[2]; mCharge[2]=nullptr;}
57  if (mXY[0]) {delete mXY[0]; mXY[0]=nullptr;}
58  if (mXY[1]) {delete mXY[1]; mXY[1]=nullptr;}
59  if (mXY[2]) {delete mXY[2]; mXY[2]=nullptr;}
60 
61  if (mEpdgeo){
62  delete mEpdgeo;
63  mEpdgeo = nullptr;
64  }
65 
66  if (mFile){
67  delete mFile;
68  mFile = nullptr;
69  }
70 
71 
72 }
73 
74 int StFcsTrackMatchMaker::Init()
75 {
76  mFcsDb = static_cast<StFcsDb *>(GetDataSet("fcsDb"));
77  if (!mFcsDb)
78  {
79  LOG_ERROR << "StFcsTrackMatchMaker::InitRun Failed to get StFcsDb" << endm;
80  return kStFatal;
81  }
82  mEpdgeo = new StEpdGeom;
83  if (mFilename)
84  {
85  LOG_INFO << "Opening " << mFilename << endm;
86  mFile = new TFile(mFilename, "RECREATE");
87 
88  mHdx[0] = new TH1F("dx_EcalTrk", "dx Ecal-Track", 100, -50, 50);
89  mHdy[0] = new TH1F("dy_EcalTrk", "dy Ecal-Track", 100, -50, 50);
90  mHdr[0] = new TH1F("dr_EcalTrk", "dr Ecal-Track", 100, 0, 20);
91  mNtrk[0] = new TH1F("NTrk_Ecal", "NTrk_Ecal", 10, 0.0, 10.0);
92  mNclu[0] = new TH1F("NEcalClu_Trk", "NEcalClu_Trk", 10, 0.0, 10.0);
93 
94  mHdx[1] = new TH1F("dx_HcalTrk", "dx Hcal-Track", 100, -50, 50);
95  mHdy[1] = new TH1F("dy_HcalTrk", "dy Hcal-Track", 100, -50, 50);
96  mHdr[1] = new TH1F("dr_HcalTrk", "dr Ecal-Track", 100, 0, 20);
97  mNtrk[1] = new TH1F("NTrk_Hcal", "NTrk_Hcal", 10, 0.0, 10.0);
98  mNclu[1] = new TH1F("NHcalClu_Trk", "NHcalClu_Trk", 10, 0.0, 10.0);
99 
100  mNtrk[2] = new TH1F("NTrk", "NTrk/evt", 50, 0.0, 100.0);
101  mNtrk[3] = new TH1F("NGoodTrk", "NGoodTrk/evt", 20, 0.0, 20.0);
102 
103  mNclu[2] = new TH1F("NEcalClu", "NEcalClu/evt", 30, 0.0, 30.0);
104  mNclu[3] = new TH1F("NHcalClu", "NHcalClu/evt", 30, 0.0, 30.0);
105 
106  mPtEt[0] = new TH1F("ETovPT_E", "ETovPT_E", 50, 0.0, 2.0);
107  mPtEt[1] = new TH1F("ETovPT_EH", "ETovPT_E+H", 50, 0.0, 2.0);
108  mPtEt2[0] = new TH2F("ETPT_E", "ETvsPT_E; ET(Ecal); TrkPT", 50, 0.0, 8.0, 50, 0.0, 8.0);
109  mPtEt2[1] = new TH2F("ETPT_EH", "ETvsPT_E+H ET(E+H); TrkPT", 50, 0.0, 8.0, 50, 0.0, 8.0);
110 
111  mCharge[0] = new TH1F("Charge_E", "Charge_E", 3, -1.5, 1.5);
112  mCharge[1] = new TH1F("Charge_H", "Charge_H", 3, -1.5, 1.5);
113  mCharge[2] = new TH1F("Charge", "Charge", 3, -1.5, 1.5);
114 
115  mXY[0] = new TH2F("XY_E", "XY_E", 50, -130, 130, 50, -110, 110);
116  mXY[1] = new TH2F("XY_H", "XY_H", 50, -130, 130, 50, -110, 110);
117  mXY[2] = new TH2F("XY", "XY", 50, -130, 130, 50, -110, 110);
118  }
119  return kStOK;
120 }
121 
123 {
124  if (mFile)
125  {
126  LOG_INFO << "Closing " << mFilename << endm;
127  mFile->Write();
128  mFile->Close();
129  }
130  return kStOK;
131 }
132 
134 {
135  StEvent *event = (StEvent *)GetInputDS("StEvent");
136  if (!event)
137  {
138  LOG_ERROR << "StFcsTrackMatchMaker::Make did not find StEvent" << endm;
139  return kStErr;
140  }
141  mFcsColl = event->fcsCollection();
142  if (!mFcsColl)
143  {
144  LOG_ERROR << "StFcsTrackMatchMaker::Make did not find StEvent->StFcsCollection" << endm;
145  return kStErr;
146  }
147  mFwdTrkColl = event->fwdTrackCollection();
148  if (!mFwdTrkColl)
149  {
150  LOG_ERROR << "StFcsTrackMatchMaker::Make did not find StEvent->fwdTrackCollection" << endm;
151  return kStErr;
152  }
153 
154  int ntrk = mFwdTrkColl->numberOfTracks();
155  int ngoodtrk = 0;
156  int nMatch[2] = {0, 0};
157  for (int itrk = 0; itrk < ntrk; itrk++)
158  {
159  StFwdTrack *trk = mFwdTrkColl->tracks()[itrk];
160  if (trk->didFitConvergeFully() == false)
161  continue;
162  ngoodtrk++;
163 
164  // get the projections
165  StFwdTrackProjection projECAL, projHCAL, projEPD;
166  LOG_INFO << "nProjections: "<< trk->mProjections.size() << endm;
167  projECAL = trk->getProjectionFor( kFcsWcalId );
168  projHCAL = trk->getProjectionFor( kFcsHcalId );
169  projEPD = trk->getProjectionFor( kFcsPresId );
170 
171  if (Debug())
172  {
173 
174  LOG_INFO << Form("Proj0: %6.2f %6.2f %6.2f", projECAL.mXYZ.x(), projECAL.mXYZ.y(), projECAL.mXYZ.z()) << endm;
175  LOG_INFO << Form("Proj1: %6.2f %6.2f %6.2f", projHCAL.mXYZ.x(), projHCAL.mXYZ.y(), projHCAL.mXYZ.z()) << endm;
176  LOG_INFO << Form("Proj2: %6.2f %6.2f %6.2f", projEPD. mXYZ.x(), projEPD. mXYZ.y(), projEPD. mXYZ.z()) << endm;
177  }
178 
179  // North or south from track
180  int ns = 0;
181  if (projECAL.mXYZ.x() > 0.0)
182  ns = 1;
183  // Look for a Ecal & Hcal match for a track
184  for (int ehp = 0; ehp < 2; ehp++)
185  {
186  StThreeVectorD proj = projECAL.mXYZ;
187  if ( ehp == 1 )
188  proj = projHCAL.mXYZ;
189 
190  int det = ehp * 2 + ns;
191  int nclu = mFcsColl->numberOfClusters(det);
192  for (int iclu = 0; iclu < nclu; iclu++)
193  {
194  StFcsCluster *clu = mFcsColl->clusters(det)[iclu];
195  float energy = clu->energy();
196  if (energy > mMinEnergy[ehp])
197  {
198  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(det, clu->x(), clu->y());
199  double dx = xyz.x() - proj.x();
200  double dy = xyz.y() - proj.y();
201  double dr = sqrt(dx * dx + dy * dy);
202  if (Debug())
203  LOG_INFO << Form("EHP=%1d dx = %6.2f - %6.2f = %6.2f dy = %6.2f - %6.2f = %6.2f dr=%6.2f",
204  ehp, xyz.x(), proj.x(), dx, xyz.y(), proj.y(), dy, dr)
205  << endm;
206  if (mFile)
207  {
208  if (fabs(dy) < mMaxDistance[ehp])
209  mHdx[ehp]->Fill(dx);
210  if (fabs(dx) < mMaxDistance[ehp])
211  mHdy[ehp]->Fill(dy);
212  mHdr[ehp]->Fill(dr);
213  }
214  if (dr < mMaxDistance[ehp])
215  {
216  if (ehp == 0)
217  {
218  trk->addEcalCluster(clu);
219  }
220  if (ehp == 1)
221  {
222  trk->addHcalCluster(clu);
223  }
224  clu->addTrack(trk);
225  nMatch[ehp]++;
226  }
227  }
228  }
229  }
230  }
231 
232  // Sort by ET/PT
233  for (int itrk = 0; itrk < ntrk; itrk++)
234  {
235  StFwdTrack *trk = mFwdTrkColl->tracks()[itrk];
236  if (trk->didFitConvergeFully() == false)
237  continue;
238  trk->sortEcalClusterByET();
239  trk->sortHcalClusterByET();
240  if (Debug())
241  {
242  LOG_INFO << Form("TRK pT=%6.2f Cg=%1d NEcal=%lu NHcal=%lu",
243  trk->momentum().perp(), trk->charge(),
244  trk->ecalClusters().size(),
245  trk->hcalClusters().size())
246  << endm;
247  }
248  }
249  for (int det = 0; det < 4; det++)
250  {
251  int ehp = det / 2;
252  int nclu = mFcsColl->numberOfClusters(det);
253  for (int iclu = 0; iclu < nclu; iclu++)
254  {
255  StFcsCluster *clu = mFcsColl->clusters(det)[iclu];
256  clu->sortTrackByPT();
257  if (Debug())
258  {
259  if (clu->energy() > mMinEnergy[ehp])
260  {
261  LOG_INFO << Form("FCS DET=%d ET=%6.2f NTrk=%ld",
262  clu->detectorId(), clu->fourMomentum().perp(), clu->tracks().size())
263  << endm;
264  }
265  }
266  }
267  }
268 
269  // Filling hitograms if file is specified
270  if (mFile)
271  {
272  mNtrk[2]->Fill(ntrk);
273  mNtrk[3]->Fill(ngoodtrk);
274  for (int itrk = 0; itrk < ntrk; itrk++)
275  {
276  StFwdTrack *trk = mFwdTrkColl->tracks()[itrk];
277  if (trk->didFitConvergeFully() == false)
278  continue;
279 
280  // get the projections
281  StFwdTrackProjection projECAL, projHCAL, projEPD;
282  projECAL = trk->getProjectionFor( kFcsWcalId );
283  projHCAL = trk->getProjectionFor( kFcsHcalId );
284  projEPD = trk->getProjectionFor( kFcsPresId );
285 
286  mCharge[2]->Fill(trk->charge());
287  mXY[2]->Fill(projECAL.mXYZ.x(), projECAL.mXYZ.y());
288  double pt = trk->momentum().perp();
289  int ne = trk->ecalClusters().size();
290  int nh = trk->hcalClusters().size();
291  mNclu[0]->Fill(double(ne));
292  mNclu[1]->Fill(double(nh));
293  double ete = 0;
294  if (ne > 0)
295  {
296  StFcsCluster *eclu = trk->ecalClusters()[0]; // Take top ET ones
297  ete = eclu->fourMomentum().perp();
298  mPtEt[0]->Fill(ete / pt);
299  mPtEt2[0]->Fill(ete, pt);
300  mCharge[0]->Fill(trk->charge());
301  mXY[0]->Fill(projECAL.mXYZ.x(), projECAL.mXYZ.y());
302  }
303  if (nh > 0)
304  {
305  StFcsCluster *hclu = trk->hcalClusters()[0]; // Take top ET ones
306  double eth = hclu->fourMomentum().perp() + ete;
307  mPtEt[1]->Fill(eth / pt);
308  mPtEt2[1]->Fill(eth, pt);
309  mCharge[1]->Fill(trk->charge());
310  mXY[1]->Fill(projHCAL.mXYZ.x(), projHCAL.mXYZ.y());
311  }
312  }
313 
314  int nc[2] = {0, 0};
315  for (int det = 0; det < 4; det++)
316  {
317  int ehp = det / 2;
318  int nclu = mFcsColl->numberOfClusters(det);
319  nc[ehp] += nclu;
320  for (int iclu = 0; iclu < nclu; iclu++)
321  {
322  StFcsCluster *clu = mFcsColl->clusters(det)[iclu];
323  if (clu->energy() > mMinEnergy[ehp])
324  {
325  mNtrk[ehp]->Fill(clu->tracks().size());
326  }
327  }
328  }
329  mNclu[2]->Fill(nc[0]);
330  mNclu[3]->Fill(nc[1]);
331 
332  LOG_INFO << Form("NTrack=%5d NGoodTrack=%3d NEcalMatch=%3d NHcalMatch=%3d", ntrk, ngoodtrk, nMatch[0], nMatch[1]) << endm;
333  }
334  return kStOK;
335 }
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
Definition: Stypes.h:40
Definition: Stypes.h:44