StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFcsPi0FinderForEcal.cxx
1 //class StFcsPi0FinderForEcal
2 //author Xilin Liang
3 //
4 //critical plots for offline QA: h1_inv_mass_cluster , h1_two_cluster_energy_allcut , h1_dgg_cluster , h2_cluster_dgg_vs_E1pE2
5 //
6 //Update 3/12/24: move the StFcsPi0FinderForEcal to StRoot , no dependent on StSpinPool , also comment out 2 unused head file
7 //
8 
9 #include "StFcsPi0FinderForEcal.h"
10 
11 #include "StEvent/StEnumerations.h"
12 #include "StEvent/StEvent.h"
13 #include "StEvent/StFcsCluster.h"
14 #include "StEvent/StFcsCollection.h"
15 #include "StEvent/StFcsHit.h"
16 #include "StEventTypes.h"
17 #include "StFcsDbMaker/StFcsDbMaker.h"
18 #include "StMessMgr.h"
19 #include "StMuDSTMaker/COMMON/StMuTypes.hh"
20 //#include "StSpinPool/StFcsQaMaker/StFcsQaMaker.h"
21 //#include "StSpinPool/StFcsRawDaqReader/StFcsRawDaqReader.h"
22 #include "StThreeVectorF.hh"
23 #include "Stypes.h"
24 #include "TBox.h"
25 #include "TCanvas.h"
26 #include "TColor.h"
27 #include "TFile.h"
28 #include "TH1.h"
29 #include "TH2.h"
30 #include "TLine.h"
31 #include "TMarker.h"
32 #include "TROOT.h"
33 #include "TString.h"
34 #include "TStyle.h"
35 #include "TText.h"
36 
37 #ifndef SKIPDefImp
38 ClassImp(StFcsPi0FinderForEcal)
39 #endif
40 
41  //------------------------
42  StFcsPi0FinderForEcal::StFcsPi0FinderForEcal(const Char_t* name) : StMaker(name) {
43 }
44 
45 StFcsPi0FinderForEcal::~StFcsPi0FinderForEcal() {}
46 
47 //-----------------------
48 Int_t StFcsPi0FinderForEcal::Init() {
49  mFcsDb = static_cast<StFcsDb*>(GetDataSet("fcsDb"));
50  if (!mFcsDb) {
51  LOG_ERROR << "StFcsEventDisplay::InitRun Failed to get StFcsDbMaker" << endm;
52  return kStFatal;
53  }
54 
55  h1_num_entries = new TH1F("h1_num_entries", "# of entries", 10, 0, 10);
56  h1_inv_mass_cluster = new TH1F("h1_inv_mass_cluster", "invariant mass plot for FCS ECal cluster", bins, m_low, m_up);
57  h1_inv_mass_cluster->SetXTitle("invariant mass [GeV]");
58  h1_Zgg_cluster = new TH1F("h1_Zgg_cluster", "Zgg", bins, 0, .7);
59  h1_Zgg_cluster->SetXTitle("Z_{gg}");
60  h1_opening_angle_cluster = new TH1F("h1_opening_angle_cluster", "FCS ECal cluster pair opening angle", bins, 0, 0.1);
61  h1_opening_angle_cluster->SetXTitle("opening angle [rad]");
62  h1_each_cluster_energy = new TH1F("h1_each_cluster_energy", "each cluster energy for FCS ECal", bins, 0, 25);
63  h1_each_cluster_energy->SetXTitle("E_{1} [GeV]");
64  h1_two_cluster_energy_nocut = new TH1F("h1_two_cluster_energy_nocut", "2 clusters energy(no cut)", bins, 0, 30);
65  h1_two_cluster_energy_nocut->SetXTitle("E_{1} + E_{2} [GeV]");
66  h1_two_cluster_energy_allcut = new TH1F("h1_two_cluster_energy_allcut", "2 clusters total energy for FCS ECal", bins, 0, 30);
67  h1_two_cluster_energy_allcut->SetXTitle("E_{1} + E_{2} [GeV]");
68  h1_dgg_cluster = new TH1F("h1_dgg_cluster", "2 cluster distance for FCS ECal", bins, 0, 250);
69  h1_dgg_cluster->SetXTitle("d_{gg} [cm]");
70  h1_Zgg_nocut_cluster = new TH1F("h1_Zgg_nocut_cluster", "Zgg without cut", bins, 0, 1);
71  h1_Zgg_nocut_cluster->SetXTitle("Z_{gg}");
72  h1_nCluster = new TH1I("h1_nCluster", "Number of clusters", 40, 0, 40);
73  h1_inv_mass_cluster_nocut = new TH1F("h1_inv_mass_cluster_nocut", "invariant mass plot (no cut)", bins, m_low, m_up);
74  h1_inv_mass_cluster_nocut->SetXTitle("invariant mass [GeV]");
75  h1_nclu_good = new TH1I("h1_nclu_good", "number of good clusters", 50, 0, 50);
76  h1_clu_nTowers = new TH1I("h1_clu_nTowers", "number of towers for cluster", 15, 0, 15);
77  h1_clu_nTowers->SetXTitle("n Towers");
78 
79  //point
80  h1_inv_mass_point = new TH1F("h1_inv_mass_point", "invariant mass plot for FCS ECal point", bins, m_low, m_up);
81  h1_inv_mass_point->SetXTitle("invariant mass [GeV]");
82  h1_Zgg_point = new TH1F("h1_Zgg_point", "energy asymmetry for FCS ECal point", bins, 0, .7);
83  h1_Zgg_point->SetXTitle("Z_{gg}");
84  h1_opening_angle_point = new TH1F("h1_opening_angle_point", "opening angle plot", bins, 0, 0.1);
85  h1_opening_angle_point->SetXTitle("opening angle [rad]");
86  h1_each_point_energy = new TH1F("h1_each_point_energy", "each point energy for FCS ECal", bins, 0, 25);
87  h1_each_point_energy->SetXTitle("E_{1} [GeV]");
88  h1_two_point_energy_nocut = new TH1F("h1_two_point_energy_nocut", "2 points energy(no cut)", bins, 0, 30);
89  h1_two_point_energy_nocut->SetXTitle("E_{1} + E_{2} [GeV]");
90  h1_two_point_energy_allcut = new TH1F("h1_two_point_energy_allcut", "2 points energy for FCS ECal point", bins, 0, 50);
91  h1_two_point_energy_allcut->SetXTitle("E_{1} + E_{2} [GeV]");
92  h1_dgg_point = new TH1F("h1_dgg_point", "d_{gg}", bins, 0, 250);
93  h1_dgg_point->SetXTitle("d_{gg} [cm]");
94  h1_Zgg_nocut_point = new TH1F("h1_Zgg_nocut_point", "Zgg without cut", bins, 0, 1);
95  h1_Zgg_nocut_point->SetXTitle("Z_{gg}");
96  h1_EcalMult_E1 = new TH1F("h1_EcalMult_E1", "Ecal Mult (E>1)", 50, 0, 50);
97  h1_EcalMult_E1->SetXTitle("Ecal Mult");
98  h1_nPoint = new TH1I("h1_nPoint", "Number of points", 40, 0, 40);
99  h1_inv_mass_point_nocut = new TH1F("h1_inv_mass_point_nocut", "#pi^{0} invariant mass plot (no cut)", bins, m_low, m_up);
100  h1_inv_mass_point_nocut->SetXTitle("#pi^{0} invariant mass [GeV]");
101  h1_npoi_good = new TH1I("h1_npoi_good", "number of good points", 50, 0, 50);
102 
103  for (int i = 0; i < 748; i++) {
104  char name_hist[50];
105  char title_hist[100];
106  sprintf(name_hist, "mass_by_Ntower_%i", i);
107  sprintf(title_hist, "invariant mass by %i North tower ", i);
108  h1list_mass_by_Ntower[i] = new TH1F(name_hist, title_hist, 50, m_low, m_up);
109  h1list_mass_by_Ntower[i]->SetXTitle("invariant mass [GeV]");
110  sprintf(name_hist, "mass_by_Stower_%i", i);
111  sprintf(title_hist, "invariant mass by %i Sorth tower ", i);
112  h1list_mass_by_Stower[i] = new TH1F(name_hist, title_hist, 50, m_low, m_up);
113  h1list_mass_by_Stower[i]->SetXTitle("invariant mass [GeV]");
114  sprintf(name_hist, "NEtower_%i", i);
115  sprintf(title_hist, "North %i tower energy spectrum", i + 1);
116  h1list_NEtower[i] = new TH1F(name_hist, title_hist, 200, 0, 20);
117  h1list_NEtower[i]->SetXTitle("Energy [GeV]");
118  sprintf(name_hist, "SEtower_%i", i);
119  sprintf(title_hist, "South %i tower energy spectrum", i + 1);
120  h1list_SEtower[i] = new TH1F(name_hist, title_hist, 200, 0, 20);
121  h1list_SEtower[i]->SetXTitle("Energy [GeV]");
122  }
123 
124  h2_EcalMult_vs_TofMult = new TH2F("h2_EcalMult_vs_TofMult", "Ecal Mult vs Tof Mult", 200, 0, 200, 200, 0, 200);
125  h2_EcalMult_vs_TofMult->SetXTitle("Tof Multiplicity");
126  h2_EcalMult_vs_TofMult->SetYTitle("Ecal Multiplicity");
127  h2_cluster_position = new TH2F("h2_cluster_position", "cluster_position", 300, -150, 150, 300, -150, 150);
128  h2_cluster_position->SetXTitle("x [cm]");
129  h2_cluster_position->SetYTitle("y [cm]");
130  h2_cluster_position_cut = new TH2F("h2_cluster_position_cut", "cluster position after cut", 400, -200, 200, 350, -150, 150);
131  h2_cluster_position_cut->SetXTitle("x [cm]");
132  h2_cluster_position_cut->SetYTitle("y [cm]");
133  h2_point_position = new TH2F("h2_point_position", "point_position", 300, -150, 150, 300, -150, 150);
134  h2_point_position->SetXTitle("x [cm]");
135  h2_point_position->SetYTitle("y [cm]");
136 
137  h2_cluster_invmass_vs_dgg = new TH2F("h2_cluster_invmass_vs_dgg", "invariant mass vs d_{gg} (cluster)", 60, 0, 60, 60, 0, 0.3);
138  h2_cluster_invmass_vs_dgg->SetXTitle("d_{gg} [cm]");
139  h2_cluster_invmass_vs_dgg->SetYTitle("m_{gg} [GeV]");
140  h2_cluster_invmass_vs_Zgg = new TH2F("h2_cluster_invmass_vs_Zgg", "invariant mass vs Z_{gg} (cluster)", 35, 0, 0.7, 60, 0, 0.3);
141  h2_cluster_invmass_vs_Zgg->SetXTitle("Z_{gg} ");
142  h2_cluster_invmass_vs_Zgg->SetYTitle("m_{gg} [GeV]");
143  h2_point_invmass_vs_dgg = new TH2F("h2_point_invmass_vs_dgg", "invariant mass vs d_{gg} (cluster)", 60, 0, 60, 60, 0, 0.6);
144  h2_point_invmass_vs_dgg->SetXTitle("d_{gg} [cm]");
145  h2_point_invmass_vs_dgg->SetYTitle("m_{gg} [GeV]");
146  h2_point_invmass_vs_Zgg = new TH2F("h2_point_invmass_vs_Zgg", "invariant mass vs Z_{gg} (cluster)", 35, 0, 0.7, 60, 0, 0.6);
147  h2_point_invmass_vs_Zgg->SetXTitle("Z_{gg} ");
148  h2_point_invmass_vs_Zgg->SetYTitle("m_{gg} [GeV]");
149  h2_cluster_dgg_vs_E1pE2 = new TH2F("h2_cluster_dgg_vs_E1pE2", "2 cluster distance vs total energy for FCS ECal cluster pair", 20, 0, 20, 60, 0, 60);
150  h2_cluster_dgg_vs_E1pE2->SetXTitle("E_{1}+E_{2} [GeV]");
151  h2_cluster_dgg_vs_E1pE2->SetYTitle("d_{gg} [cm]");
152  h2_point_dgg_vs_E1pE2 = new TH2F("h2_point_dgg_vs_E1pE2", "energy asymmetry vs total energy for FCS ECal point pair", 20, 0, 20, 60, 0, 60);
153  h2_point_dgg_vs_E1pE2->SetXTitle("E_{1}+E_{2} [GeV]");
154  h2_point_dgg_vs_E1pE2->SetYTitle("d_{gg} [cm]");
155 
156  return kStOK;
157 }
158 
159 //-----------------------
161  if (filename.length() == 0) return kStOk;
162  const char* fn = filename.c_str();
163  TFile MyFile(fn, "RECREATE");
164  h1_num_entries->Write();
165  h1_inv_mass_cluster->Write();
166  h1_Zgg_cluster->Write();
167  h1_opening_angle_cluster->Write();
168  h1_each_cluster_energy->Write();
169  h1_two_cluster_energy_nocut->Write();
170  h1_two_cluster_energy_allcut->Write();
171  h1_dgg_cluster->Write();
172  h1_Zgg_nocut_cluster->Write();
173  h1_inv_mass_cluster_nocut->Write();
174  h1_nCluster->Write();
175  h1_nclu_good->Write();
176  h1_clu_nTowers->Write();
177 
178  h1_num_entries->Write();
179  h1_inv_mass_point->Write();
180  h1_Zgg_point->Write();
181  h1_opening_angle_point->Write();
182  h1_each_point_energy->Write();
183  h1_two_point_energy_nocut->Write();
184  h1_two_point_energy_allcut->Write();
185  h1_dgg_point->Write();
186  h1_Zgg_nocut_point->Write();
187  h1_EcalMult_E1->Write();
188  h1_inv_mass_point_nocut->Write();
189  h1_nPoint->Write();
190  h1_npoi_good->Write();
191 
192  for (int i = 0; i < 748; i++) {
193  h1list_mass_by_Ntower[i]->Write();
194  h1list_mass_by_Stower[i]->Write();
195  h1list_NEtower[i]->Write();
196  h1list_SEtower[i]->Write();
197  }
198 
199  h2_EcalMult_vs_TofMult->Write();
200  h2_cluster_position->Write();
201  h2_cluster_position_cut->Write();
202  h2_point_position->Write();
203  h2_cluster_invmass_vs_dgg->Write();
204  h2_cluster_invmass_vs_Zgg->Write();
205  h2_point_invmass_vs_dgg->Write();
206  h2_point_invmass_vs_Zgg->Write();
207  h2_cluster_dgg_vs_E1pE2->Write();
208  h2_point_dgg_vs_E1pE2->Write();
209 
210  MyFile.Close();
211  return kStOK;
212 }
213 
214 //----------------------
216  cout << "Start using StFcsPi0FinderForECal" << endl;
217  StEvent* event = (StEvent*)GetInputDS("StEvent");
218  if (!event) {
219  LOG_ERROR << "StFcsPi0FinderForEcal::Make did not find StEvent" << endm;
220  return kStErr;
221  }
222  mFcsColl = event->fcsCollection();
223  if (!mFcsColl) {
224  LOG_ERROR << "StFcsPi0FinderForEcal::Make did not find StEvent->StFcsCollection" << endm;
225  return kStErr;
226  }
227 
228  if (mNAccepted < mMaxEvents) {
229  mNEvents++;
230  cout << "current event:" << mNEvents << endl;
231  if (mFilter == 1 && mFcsColl->numberOfHits(0) + mFcsColl->numberOfHits(1) + mFcsColl->numberOfHits(2) + mFcsColl->numberOfHits(3) == 0) return kStOK;
232 
233  //TOF mult cut
234  int tofMult = 0;
235  const StTriggerData* trgdata = event->triggerData();
236  if(!trgdata && StMuDst::event()) trgdata = StMuDst::event()->triggerData();
237  if(trgdata){
238  tofMult = trgdata->tofMultiplicity();
239  LOG_DEBUG<<"TOF mult="<<tofMult<<endm;
240  if (tofMult > 100) return kStOK;
241  }else{
242  LOG_WARN << "No TriggerData found in StEvent nor Mudst. No TOFMult cut"<<endm;
243  }
244 
245  //TPC ZVERTEX
246  float zTPC=-999.0;
247  StPrimaryVertex* tpcvtx = event->primaryVertex();
248  if(tpcvtx) {
249  zTPC=tpcvtx->position().z();
250  }else{
251  if (StMuDst::numberOfPrimaryVertices() > 0){
253  if(mutpcvtx) zTPC=mutpcvtx->position().z();
254  }
255  }
256 
257  //BBC ZVERTEX
258  float zBBC=-999.0;
259  if(trgdata) zBBC = (4096 - trgdata->bbcTimeDifference())*0.016*30.0/2.0;
260  if(zBBC<-200 || zBBC>200) zBBC=-999;
261 
262  //VPD ZVERTEX from MuDst(TOF data)
263  float zVPD=-999.0;
264  if(StMuDst::btofHeader()) zVPD=StMuDst::btofHeader()->vpdVz();
265 
266  LOG_INFO << Form("ZTPX = %6.2f ZBBC = %6.2f ZVPD = %6.2f",zTPC,zBBC,zVPD) << endm;
267 
268  //test getLorentzVector
269  StThreeVectorD xyz(20,0,720);
270  StLorentzVectorD pbbc,ptpc,pvpd;
271  StLorentzVectorD p0 = mFcsDb->getLorentzVector((xyz), 10, 0); LOG_INFO << "Zero " << p0 << endm;
272  if(zBBC>-200) {pbbc = mFcsDb->getLorentzVector((xyz), 10, zBBC); LOG_INFO << "BBC " << pbbc << endm;}
273  if(zTPC>-200) {ptpc = mFcsDb->getLorentzVector((xyz), 10, zTPC); LOG_INFO << "TPC " << ptpc << endm;}
274  if(zVPD>-200) {pvpd = mFcsDb->getLorentzVector((xyz), 10, zVPD); LOG_INFO << "VPD " << pvpd << endm;}
275 
276  mNAccepted++;
277  int total_nc = 0;
278  int total_np = 0;
279  int check_fillclu = 0;
280  int check_fillpnt = 0;
281 
282  int n_EcalMult = 0;
283  int n_EcalClustMult = 0;
284 
285  int n_EcalClust_cut = 0;
286  int n_EcalPoint_cut = 0;
287  int n_Ecal_cut = 0;
288 
289  float bestclu_invmass = -1;
290  float bestclu_totalE = -1;
291  float bestclu_dgg = -1;
292  float bestclu_Zgg = -1;
293  float bestclu_opening_angle = -1;
294  float bestpnt_invmass = -1;
295  float bestpnt_totalE = -1;
296  float bestpnt_dgg = -1;
297  float bestpnt_Zgg = -1;
298  float bestpnt_opening_angle = -1;
299  int best_tower_id1_cluster = -1;
300  int best_tower_id2_cluster = -1;
301  int best_tower_det_cluster = -1;
302 
303  for (int det = 0; det < 2; det++) {
304  int check_Ecal = mFcsDb->ecalHcalPres(det);
305  if (check_Ecal != 0) continue;
306 
307  StSPtrVecFcsCluster& clusters = mFcsColl->clusters(det);
308  int nc = mFcsColl->numberOfClusters(det);
309 
310  total_nc = nc + total_nc;
311 
312  StSPtrVecFcsPoint& points = mFcsColl->points(det);
313  int np = mFcsColl->numberOfPoints(det);
314 
315  total_np = np + total_np;
316 
317  StSPtrVecFcsHit& hits = mFcsColl->hits(det);
318  int nh = mFcsColl->numberOfHits(det);
319  if (mDebug > 0) LOG_INFO << Form("StFcsEventDisplay Det=%1d nhit=%4d nclu=%3d", det, nh, nc) << endm;
320 
321  for (int i = 0; i < nh; i++) {
322  StFcsHit* hit = hits[i];
323  unsigned short hit_id = hit->id();
324 
325  float hit_energy = hit->energy();
326  if (det == 0) {
327  h1list_NEtower[hit_id]->Fill(hit_energy);
328  } else if (det == 1) {
329  h1list_SEtower[hit_id]->Fill(hit_energy);
330  } //fill in energy spectrum for tower
331  if (hit_energy > E_min) n_Ecal_cut++;
332  n_EcalMult++;
333  }
334 
335  // //no cut (cluster)
336  for (int i = 0; i < nc; i++) {
337  StFcsCluster* clu = clusters[i];
338  float clu_energy = clu->energy();
339  float clu_x = clu->x();
340  float clu_y = clu->y();
341  StThreeVectorD cluPos = mFcsDb->getStarXYZfromColumnRow(det, clu_x, clu_y);
342  float cluPos_x = cluPos.x();
343  float cluPos_y = cluPos.y();
344 
345  n_EcalClustMult++;
346  if ((clu_energy > E_min)) {
347  n_EcalClust_cut++;
348  }
349 
350  h2_cluster_position->Fill(cluPos_x, cluPos_y);
351  h1_each_cluster_energy->Fill(clu_energy);
352  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(det, clu_x, clu_y);
353  StLorentzVectorD p = mFcsDb->getLorentzVector((xyz), clu_energy, 0);
354  if (i == nc - 1) continue;
355  for (int j = i + 1; j < nc; j++) {
356  StFcsCluster* cluj = clusters[j];
357  float cluj_energy = cluj->energy();
358  float cluj_x = cluj->x();
359  float cluj_y = cluj->y();
360  StThreeVectorD clujPos = mFcsDb->getStarXYZfromColumnRow(det, cluj_x, cluj_y);
361 
362  h1_two_cluster_energy_nocut->Fill(clu_energy + cluj_energy);
363  float zgg = (abs(clu_energy - cluj_energy)) / (clu_energy + cluj_energy);
364  h1_Zgg_nocut_cluster->Fill(zgg);
365  StThreeVectorD xyzj = mFcsDb->getStarXYZfromColumnRow(det, cluj->x(), cluj->y());
366  StLorentzVectorD pj = mFcsDb->getLorentzVector((xyzj), cluj->energy(), 0);
367  h1_inv_mass_cluster_nocut->Fill((p + pj).m());
368  }
369  }
370 
371  // //no cut (point)
372  for (int i = 0; i < np; i++) {
373  StFcsPoint* pnt = points[i];
374  float pnt_x = pnt->x();
375  float pnt_y = pnt->y();
376  float pnt_energy = pnt->energy();
377  if ((pnt_energy > E_min)) {
378  n_EcalPoint_cut++;
379  }
380  StThreeVectorD poiPos = mFcsDb->getStarXYZfromColumnRow(det, pnt_x, pnt_y);
381  float poiPos_x = poiPos.x();
382  float poiPos_y = poiPos.y();
383  h2_point_position->Fill(poiPos_x, poiPos_y);
384  h1_each_point_energy->Fill(pnt_energy);
385  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(det, pnt_x, pnt_y);
386  StLorentzVectorD p = mFcsDb->getLorentzVector((xyz), pnt_energy, 0);
387  if (i == np - 1) continue;
388  for (int j = i + 1; j < np; j++) {
389  StFcsPoint* pntj = points[j];
390  float pntj_energy = pntj->energy();
391 
392  h1_two_point_energy_nocut->Fill(pnt_energy + pntj_energy);
393  float zgg = (abs(pnt_energy - pntj_energy)) / (pnt_energy + pntj_energy);
394  h1_Zgg_nocut_point->Fill(zgg);
395  StThreeVectorD xyzj = mFcsDb->getStarXYZfromColumnRow(det, pntj->x(), pntj->y());
396  StLorentzVectorD pj = mFcsDb->getLorentzVector((xyzj), pntj->energy(), 0);
397  h1_inv_mass_point_nocut->Fill((p + pj).m());
398  }
399  }
400 
401  //start cut (cluster)
402  for (int i = 0; i < nc; i++) {
403  StFcsCluster* clu = clusters[i];
404  float clu_x = clu->x();
405  float clu_y = clu->y();
406  float clu_energy = clu->energy();
407  h1_clu_nTowers->Fill(clu->nTowers());
408  if ((clu->energy()) < E_min) continue;
409 
410  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(det, clu_x, clu_y);
411  StLorentzVectorD p = mFcsDb->getLorentzVector((xyz), clu_energy, 0);
412  h2_cluster_position_cut->Fill(xyz.x(), xyz.y());
413  for (int j = i + 1; j < nc; j++) {
414  StFcsCluster* cluj = clusters[j];
415  float cluj_energy = cluj->energy();
416 
417  if ((cluj->energy()) < E_min) continue;
418  float zgg = (fabs(clu_energy - cluj_energy)) / (clu_energy + cluj_energy);
419  if (zgg > 0.7) continue;
420  float cluj_x = cluj->x();
421  float cluj_y = cluj->y();
422  StThreeVectorD xyzj = mFcsDb->getStarXYZfromColumnRow(det, cluj_x, cluj_y);
423  StLorentzVectorD pj = mFcsDb->getLorentzVector((xyzj), cluj_energy, 0);
424 
425  if ((clu_energy + cluj_energy) > bestclu_totalE) {
426  check_fillclu = 1;
427  bestclu_invmass = ((p + pj).m());
428  bestclu_totalE = (clu_energy + cluj_energy);
429  bestclu_dgg = (sqrt((xyz[0] - xyzj[0]) * (xyz[0] - xyzj[0]) + (xyz[1] - xyzj[1]) * (xyz[1] - xyzj[1]) + (xyz[2] - xyzj[2]) * (xyz[2] - xyzj[2])));
430  bestclu_Zgg = fabs((clu_energy - cluj_energy) / (clu_energy + cluj_energy));
431  bestclu_opening_angle = acos((xyz[0] * xyzj[0] + xyz[1] * xyzj[1] + xyz[2] * xyzj[2]) / (sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1] + xyz[2] * xyz[2]) * sqrt(xyzj[0] * xyzj[0] + xyzj[1] * xyzj[1] + xyzj[2] * xyzj[2])));
432 
433  int pnti_nTowers = clu->nTowers();
434  int pntj_nTowers = cluj->nTowers();
435  StPtrVecFcsHit& clui_hits = clu->hits();
436  StPtrVecFcsHit& cluj_hits = cluj->hits();
437  float max_energy_tower = -1;
438  unsigned short tower_id = 0;
439  for (int k = 0; k < pnti_nTowers; k++) {
440  if (clui_hits[k]->energy() > max_energy_tower) {
441  max_energy_tower = clui_hits[k]->energy();
442  tower_id = clui_hits[k]->id();
443  best_tower_det_cluster = det;
444  }
445  }
446  best_tower_id1_cluster = tower_id;
447  max_energy_tower = -1;
448  for (int k = 0; k < pntj_nTowers; k++) {
449  if ((cluj_hits[k]->energy()) > max_energy_tower) {
450  max_energy_tower = cluj_hits[k]->energy();
451  tower_id = cluj_hits[k]->id();
452  }
453  }
454  best_tower_id2_cluster = tower_id;
455  }
456  }
457  }
458 
459  //
460  //
461  // //all cut (point)
462  for (int i = 0; i < np - 1; i++) {
463  StFcsPoint* pnt = points[i];
464  float pnt_x = pnt->x();
465  float pnt_y = pnt->y();
466  float pnt_energy = pnt->energy();
467  if (pnt_energy < E_min) continue;
468  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(det, pnt_x, pnt_y);
469  StLorentzVectorD p = mFcsDb->getLorentzVector((xyz), pnt_energy, 0);
470 
471  for (int j = i + 1; j < np; j++) {
472  StFcsPoint* pntj = points[j];
473  float pntj_energy = pntj->energy();
474 
475  if ((pntj->energy()) < E_min) continue;
476  float zgg = (abs(pnt_energy - pntj_energy)) / (pnt_energy + pntj_energy);
477  if (zgg >= 0.7) continue;
478  StThreeVectorD xyzj = mFcsDb->getStarXYZfromColumnRow(det, pntj->x(), pntj->y());
479  StLorentzVectorD pj = mFcsDb->getLorentzVector((xyzj), pntj->energy(), 0);
480  if ((pnt_energy + pntj_energy) > bestpnt_totalE) {
481  check_fillpnt = 1;
482  bestpnt_invmass = ((p + pj).m());
483  bestpnt_totalE = (pnt_energy + pntj_energy);
484  bestpnt_dgg = (sqrt((xyz[0] - xyzj[0]) * (xyz[0] - xyzj[0]) + (xyz[1] - xyzj[1]) * (xyz[1] - xyzj[1]) + (xyz[2] - xyzj[2]) * (xyz[2] - xyzj[2])));
485  bestpnt_Zgg = fabs((pnt_energy - pntj_energy) / (pnt_energy + pntj_energy));
486  bestpnt_opening_angle = acos((xyz[0] * xyzj[0] + xyz[1] * xyzj[1] + xyz[2] * xyzj[2]) / (sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1] + xyz[2] * xyz[2]) * sqrt(xyzj[0] * xyzj[0] + xyzj[1] * xyzj[1] + xyzj[2] * xyzj[2])));
487  }
488  }
489  }
490  }
491 
492  h1_nCluster->Fill(total_nc);
493  h1_nclu_good->Fill(n_EcalClust_cut);
494  h1_nPoint->Fill(total_np);
495  h1_npoi_good->Fill(n_EcalPoint_cut);
496  h2_EcalMult_vs_TofMult->Fill(tofMult, n_Ecal_cut);
497  if (n_Ecal_cut > 40) {
498  return kStOK;
499  }
500  if (check_fillclu == 1) {
501  h1_inv_mass_cluster->Fill(bestclu_invmass);
502  h1_two_cluster_energy_allcut->Fill(bestclu_totalE);
503  h2_cluster_invmass_vs_dgg->Fill(bestclu_dgg, bestclu_invmass);
504  h2_cluster_invmass_vs_Zgg->Fill(bestclu_Zgg, bestclu_invmass);
505  h1_Zgg_cluster->Fill(bestclu_Zgg);
506  h1_opening_angle_cluster->Fill(bestclu_opening_angle);
507  h1_dgg_cluster->Fill(bestclu_dgg);
508  h2_cluster_dgg_vs_E1pE2->Fill(bestclu_totalE, bestclu_dgg);
509  if (best_tower_det_cluster == 0) {
510  h1list_mass_by_Ntower[best_tower_id1_cluster]->Fill(bestclu_invmass);
511  h1list_mass_by_Ntower[best_tower_id2_cluster]->Fill(bestclu_invmass);
512  }
513  if (best_tower_det_cluster == 1) {
514  h1list_mass_by_Stower[best_tower_id1_cluster]->Fill(bestclu_invmass);
515  h1list_mass_by_Stower[best_tower_id2_cluster]->Fill(bestclu_invmass);
516  }
517  }
518 
519  if (check_fillpnt == 1) {
520  h1_inv_mass_point->Fill(bestpnt_invmass);
521  h1_two_point_energy_allcut->Fill(bestpnt_totalE);
522  h2_point_invmass_vs_dgg->Fill(bestpnt_dgg, bestpnt_invmass);
523  h2_point_invmass_vs_Zgg->Fill(bestpnt_Zgg, bestpnt_invmass);
524  h1_Zgg_point->Fill(bestpnt_Zgg);
525  h1_opening_angle_point->Fill(bestpnt_opening_angle);
526  h1_dgg_point->Fill(bestpnt_dgg);
527  h2_point_dgg_vs_E1pE2->Fill(bestpnt_totalE, bestpnt_dgg);
528  }
529  }
530  return kStOK;
531 }
static StMuPrimaryVertex * primaryVertex()
return pointer to current primary vertex
Definition: StMuDst.h:322
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
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
int ecalHcalPres(int det) const
Ecal North=0, Ecal South=1, Hcal North=2, Hcal South=3, Pres=4/5.
Definition: StFcsDb.cxx:423
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
Definition: Stypes.h:40
static StBTofHeader * btofHeader()
returns pointer to the btofHeader - dongx
Definition: StMuDst.h:423
Definition: Stypes.h:44
Definition: Stypes.h:41