StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFcsPi0ReconstructionMaker.cxx
1 //class StFcsPi0ReconstructionMaker
2 //author Xilin Liang
3 //
4 //April 2, 2020:
5 //
6 
7 #include "StFcsPi0ReconstructionMaker.h"
8 #include "StEvent/StEnumerations.h"
9 #include "StMessMgr.h"
10 #include "Stypes.h"
11 #include "StMuDSTMaker/COMMON/StMuTypes.hh"
12 
13 #include "StThreeVectorF.hh"
14 #include "StFcsDbMaker/StFcsDb.h"
15 #include "StRoot/StSpinPool/StFcsQaMaker/StFcsQaMaker.h"
16 #include "StRoot/StSpinPool/StFcsRawDaqReader/StFcsRawDaqReader.h"
17 #include "StRoot/StEvent/StTriggerData.h"
18 #include "StRoot/StEvent/StFcsCollection.h"
19 #include "StRoot/StEvent/StFcsHit.h"
20 #include "StRoot/StEvent/StFcsCluster.h"
21 #include "StRoot/StEvent/StFcsPoint.h"
22 #include "StRoot/StEvent/StEvent.h"
23 #include "StEventTypes.h"
24 #include "TCanvas.h"
25 #include "TFile.h"
26 #include "TH1F.h"
27 #include "TH2F.h"
28 #include "TLine.h"
29 #include "TBox.h"
30 #include "TMarker.h"
31 #include "TString.h"
32 #include "TColor.h"
33 #include "TStyle.h"
34 #include "TText.h"
35 #include "TROOT.h"
36 
37 ClassImp(StFcsPi0ReconstructionMaker);
38 
39 //------------------------
40 StFcsPi0ReconstructionMaker::StFcsPi0ReconstructionMaker(const Char_t* name):
41  StMaker(name){}
42 
43 StFcsPi0ReconstructionMaker::~StFcsPi0ReconstructionMaker(){}
44 
45 //-----------------------
46 Int_t StFcsPi0ReconstructionMaker::Init()
47 {
48  mFcsDb=static_cast<StFcsDb*>(GetDataSet("fcsDb"));
49  if(!mFcsDb){
50  LOG_ERROR << "StFcsEventDisplay::InitRun Failed to get StFcsDbMaker" << endm;
51  return kStFatal;
52  }
53 
54  h0=new TH1F("h0","# of entries",10,0,10); //check the number of event;
55 
56  h2_EcalMult_vs_TOFMult=new TH2F("h2_EcalMult_vs_TOFMult","Ecal Multiplicity vs TOF Multiplicity",200,-.5,199.5,71,-.5,70.5);//Ecal Multiplicity vs TOF Multiplicity
57  h2_EcalMult_vs_TOFMult->SetXTitle("tof Multiplicity");
58  h2_EcalMult_vs_TOFMult->SetYTitle("Ecal Multiplicity");
59 
60  h2_EcalMult_vs_TOFMult_E_1=new TH2F("h2_EcalMult_vs_TOFMult_E_1","Ecal Multiplicity vs TOF Multiplicity",200,-.5,199.5,70,-.5,69.5);//Ecal Multiplicity vs TOF Multiplicity after counting tower energy > 1 GeV
61  h2_EcalMult_vs_TOFMult_E_1->SetXTitle("TOF Multiplicity");
62  h2_EcalMult_vs_TOFMult_E_1->SetYTitle("Ecal Multiplicity(E>1)");
63 
64  //for cluster finder
65  h1_inv_mass_cluster=new TH1F("h1_inv_mass_cluster","invariant mass plot (for cluster finder analysis)",100,m_low,m_up); // invariant mass for cluster finder analysis
66  h1_inv_mass_cluster->SetXTitle("invariant mass [GeV]");
67 
68  h1_Zgg_cluster=new TH1F("h1_Zgg_cluster","Zgg (cluster)",70,0,.7); // energy asymmetry for cluster finder analysis, define as : (abs(E1 - E2) / (E1 + E2)
69  h1_Zgg_cluster->SetXTitle("Z_{gg}");
70 
71  h1_opening_angle_cluster=new TH1F("h1_opening_angle_cluster","opening angle plot (for cluster finder analysis)",bins,0,0.1); // opening angle
72  h1_opening_angle_cluster->SetXTitle("opening angle [rad]");
73 
74  h1_each_cluster_energy_nocut=new TH1F("h1_each_cluster_energy_nocut","each cluster energy(without cut for cluster finder)",200,0,200); // store each cluster energy (without cut for cluster finder)
75  h1_each_cluster_energy_nocut->SetXTitle("E_{1} [GeV]");
76 
77  h1_each_cluster_energy_aftercut=new TH1F("h1_each_cluster_energy_aftercut","each cluster energy(after cut for cluster finder)",200,0,200);//each cluster energy after cut
78  h1_each_cluster_energy_aftercut->SetXTitle("E_{1} [GeV]");
79 
80  h1_two_cluster_energy_nocut=new TH1F("h1_two_cluster_energy_nocut","E_{1} + E_{2} (without cut for cluster finder)",300,0,300); //E1 +E2 without cut
81  h1_two_cluster_energy_nocut->SetXTitle("E_{1} + E_{2} [GeV]");
82 
83  h1_two_cluster_energy_aftercut=new TH1F("h1_two_cluster_energy_aftercut","E_{1} + E_{2} (all cuts for cluster finder)",300,2,300);//E1 + E2 after all cuts
84  h1_two_cluster_energy_aftercut->SetXTitle("E_{1} + E_{2} [GeV]");
85 
86  h1_dgg_cluster=new TH1F("h1_dgg_cluster","d_{gg} (cluster)",80,0,60);//distance bewteen two clusters at the Ecal
87  h1_dgg_cluster->SetXTitle("d_{gg} [cm]");
88 
89  h2_EcalClustMult_vs_TOFMult=new TH2F("h2_EcalClustMult_vs_TOFMult","Ecal cluster Multiplicity vs tof Multiplicity",200,-.5,199.5,40,-.5,39.5);
90  h2_EcalClustMult_vs_TOFMult->SetXTitle("tof Multiplicity");
91  h2_EcalClustMult_vs_TOFMult->SetYTitle("Ecal Cluster Multiplicity");
92 
93  h2_EcalClustMult_vs_TOFMult_E_1=new TH2F("h2_EcalClustMult_vs_TOFMult_E_1","Ecal cluster Multiplicity vs TOF Multiplicity",50,-.5,49.5,40,-.5,39.5); //Ecal cluster Multiplicity vs TOF Multiplicity after counting cluster energy > 1 GeV
94  h2_EcalClustMult_vs_TOFMult_E_1->SetXTitle("TOF Multiplicity");
95  h2_EcalClustMult_vs_TOFMult_E_1->SetYTitle("Ecal Cluster Multiplicity(E>1)");
96 
97  h2_dgg_cluster_vs_E1pE2_cluster=new TH2F("h2_dgg_cluster_vs_E1pE2_cluster","d_{gg} vs E_{1} + E_{2} for cluster finder",45,8,25,60,-.5,59.5);
98  h2_dgg_cluster_vs_E1pE2_cluster->SetYTitle("d_{gg}");
99  h2_dgg_cluster_vs_E1pE2_cluster->SetXTitle("E_{1} + E_{2} [GeV]");
100 
101  h2_mass_cluster_vs_Zgg_cluster=new TH2F("h2_mass_cluster_vs_Zgg_cluster","invariant mass vs Z_{gg} for cluster finder",70,0,.7,80,0,.8);
102  h2_mass_cluster_vs_Zgg_cluster->SetXTitle("Z_{gg}");
103  h2_mass_cluster_vs_Zgg_cluster->SetYTitle("invariant mass [GeV]");
104 
105  h2_cluster_position=new TH2F("h2_cluster_position","cluster position y vs x (cluster finder)", 80,55,135,80,-105,-25);
106  h2_cluster_position->SetXTitle("x [cm]");
107  h2_cluster_position->SetYTitle("y [cm]");
108 
109  h2_mass_cluster_vs_dgg_cluster=new TH2F("h2_mass_cluster_vs_dgg_cluster","invariant mass vs d_{gg} (cluster finder)",45,5,50,50,m_low,m_up);
110  h2_mass_cluster_vs_dgg_cluster->SetXTitle("d_{gg} [cm]");
111  h2_mass_cluster_vs_dgg_cluster->SetYTitle("invariant mass [GeV]");
112 
113  //for point maker
114  h1_inv_mass_point=new TH1F("h1_inv_mass_point","invariant mass plot for point maker",100,m_low,m_up);//invariant mass for point maker
115  h1_inv_mass_point->SetXTitle("invariant mass [GeV]");
116 
117  h1_Zgg_point=new TH1F("h1_Zgg_point","Zgg (point)",70,0,.7);//energy asymmetry for point maker analysis, define as : (abs(E1 - E2) / (E1 + E2)
118  h1_Zgg_point->SetXTitle("Z_{gg}");
119 
120  h1_opening_angle_point=new TH1F("h1_opening_angle_point","opening angle plot (point)",bins,0,0.1);//opening angle for point maker
121  h1_opening_angle_point->SetXTitle("opening angle [rad]");
122 
123  h1_each_point_energy_nocut=new TH1F("h1_each_point_energy_nocut","each point energy(no cut) (point)",200,0,200);//each point energy without cut
124  h1_each_point_energy_nocut->SetXTitle("E_{1} [GeV]");
125 
126  h1_each_point_energy_aftercut=new TH1F("h1_each_point_energy_aftercut","each point energy(after cut) (point)",200,0,200);//each point energy after cuts
127  h1_each_point_energy_aftercut->SetXTitle("E_{1} [GeV]");
128 
129  h1_two_point_energy_nocut=new TH1F("h1_two_point_energy_nocut","E_{1} + E_{2} without cut for point maker",300,0,300);
130  h1_two_point_energy_nocut->SetXTitle("E_{1} + E_{2} [GeV]");
131 
132  h1_two_point_energy_aftercut=new TH1F("h1_two_point_energy_aftercut","E_{1} + E_{2} after cuts for point maker",300,2,300);
133  h1_two_point_energy_aftercut->SetXTitle("E_{1} + E_{2} [GeV]");
134 
135  h1_dgg_point=new TH1F("h1_dgg_point","d_{gg} (point)",80,0,60);//distance bewteen two points at Ecal.
136  h1_dgg_point->SetXTitle("d_{gg} [cm]");
137 
138  h1_x_rPosition_point=new TH1F("h1_x_rPosition_point","point position(x) in cell",50,0.0,1.0);//point position(x) in cell
139  h1_x_rPosition_point->SetXTitle("point position (x in rPosition) [cm]");
140  h1_y_rPosition_point=new TH1F("h1_y_rPosition_point","point position(y) in cell",50,0.0,1.0);//point position(y) in cell
141  h1_y_rPosition_point->SetXTitle("point position (y in rPosition) [cm]");
142 
143  h2_EcalPointMult_vs_TOFMult=new TH2F("h2_EcalPointMult_vs_TOFMult","Ecal point Multiplicity vs tof Multiplicity",200,-.5,199.5,40,-.5,39.5);//Ecal Point Multiplicity us TOF Multiplicity
144  h2_EcalPointMult_vs_TOFMult->SetXTitle("tof Multiplicity");
145  h2_EcalPointMult_vs_TOFMult->SetYTitle("Ecal Pointer Multiplicity");
146 
147  h2_EcalPointMult_vs_TOFMult_E_1=new TH2F("h2_EcalPointMult_vs_TOFMult_E_1","Ecal point Multiplicity vs TOF Multiplicity",200,-.5,199.5,40,-.5,39.5);//Ecal point Multiplicity vs TOF Multiplicity after counting point energy > 1 GeV
148  h2_EcalPointMult_vs_TOFMult_E_1->SetXTitle("TOF Multiplicity");
149  h2_EcalPointMult_vs_TOFMult_E_1->SetYTitle("Ecal Pointer Multiplicity(E>1)");
150 
151  h2_dgg_point_vs_E1pE2_point=new TH2F("h2_dgg_point_vs_E1pE2_point","d_{gg} vs E_{1} + E_{2} for point maker",45,8,20,60,-.5,59.5);
152  h2_dgg_point_vs_E1pE2_point->SetYTitle("d_{gg}");
153  h2_dgg_point_vs_E1pE2_point->SetXTitle("E_{1} + E_{2} [GeV]");
154  h2_mass_point_vs_Zgg_point=new TH2F("h2_mass_point_vs_Zgg_point","invariant mass vs Z_{gg} for point maker",70,0,.7,80,0,.8);
155  h2_mass_point_vs_Zgg_point->SetXTitle("Z_{gg}");
156  h2_mass_point_vs_Zgg_point->SetYTitle("invariant mass [GeV]");
157  h2_point_position=new TH2F("h2_point_position","point position y vs x for point maker", 80,55,135,80,-105,-25);
158  h2_point_position->SetXTitle("x [cm]");
159  h2_point_position->SetYTitle("y [cm]");
160 
161  h2_mass_point_vs_dgg_point=new TH2F("h2_mass_point_vs_dgg_point","invariant mass vs d_{gg} for point maker",45,5,50,50,m_low,m_up);
162  h2_mass_point_vs_dgg_point->SetXTitle("d_{gg} [cm]");
163  h2_mass_point_vs_dgg_point->SetYTitle("invariant mass [GeV]");
164 
165  //for towers
166  for (int i=0;i<64;i++)
167  {
168  char name_hist[50];
169  char title_hist[100];
170  sprintf(name_hist,"mass_cluster_by_tower_%i",i);
171  sprintf(title_hist,"invariant mass sorted by %i tower for cluster finder",i);
172  h1list_mass_cluster_by_tower[i]=new TH1F(name_hist,title_hist,40,m_low,m_up);//invariant mass sorted by #i tower for cluster finder
173  h1list_mass_cluster_by_tower[i]->SetXTitle("invariant mass [GeV]");
174 
175  sprintf(name_hist,"mass_point_by_tower_%i",i);
176  sprintf(title_hist,"invariant mass sorted by %i tower (point)",i);
177  h1list_mass_point_by_tower[i]=new TH1F(name_hist,title_hist,40,m_low,m_up);//invariant mass sorted by #i tower for point maker
178  h1list_mass_point_by_tower[i]->SetXTitle("invariant mass [GeV]");
179 
180  sprintf(name_hist,"Etower_%i",i);
181  sprintf(title_hist,"%i tower energy spectrum",i);
182  h1list_Etower[i]=new TH1F(name_hist, title_hist,200,0,200);//energy spectrum for towers
183  h1list_Etower[i]->SetXTitle("Energy [GeV]");
184 
185  }
186 
187  return kStOK;
188 
189  }
190 
191 
192  //-----------------------
193  Int_t StFcsPi0ReconstructionMaker::Finish()
194  {
195  const char* fn=filename.c_str();
196  TFile* MyFile=new TFile(fn,"RECREATE");
197  h0->Write();
198  h1_inv_mass_cluster->Write();
199  h1_Zgg_cluster->Write();
200  h1_opening_angle_cluster->Write();
201  h1_each_cluster_energy_nocut->Write();
202  h1_each_cluster_energy_aftercut->Write();
203  h1_two_cluster_energy_nocut->Write();
204  h1_two_cluster_energy_aftercut->Write();
205  h1_dgg_cluster->Write();
206  h2_EcalMult_vs_TOFMult->Write();
207  h2_EcalMult_vs_TOFMult_E_1->Write();
208  h2_EcalClustMult_vs_TOFMult->Write();
209  h2_EcalClustMult_vs_TOFMult_E_1->Write();
210  h2_dgg_cluster_vs_E1pE2_cluster->Write();
211  h2_mass_cluster_vs_Zgg_cluster->Write();
212  h2_cluster_position->Write();
213  h2_mass_cluster_vs_dgg_cluster->Write();
214  h1_inv_mass_point->Write();
215  h1_Zgg_point->Write();
216  h1_opening_angle_point->Write();
217  h1_each_point_energy_nocut->Write();
218  h1_each_point_energy_aftercut->Write();
219  h1_two_point_energy_nocut->Write();
220  h1_two_point_energy_aftercut->Write();
221  h1_dgg_point->Write();
222  h1_x_rPosition_point->Write();
223  h1_y_rPosition_point->Write();
224  h2_EcalPointMult_vs_TOFMult->Write();
225  h2_EcalPointMult_vs_TOFMult_E_1->Write();
226  h2_dgg_point_vs_E1pE2_point->Write();
227  h2_mass_point_vs_Zgg_point->Write();
228  h2_point_position->Write();
229  h2_mass_point_vs_dgg_point->Write();
230  for (int i=0;i<64;i++)
231  {
232  h1list_mass_cluster_by_tower[i]->Write();
233  h1list_mass_point_by_tower[i]->Write();
234  h1list_Etower[i]->Write();
235  }
236  MyFile->Close();
237  return kStOK;
238  }
239 
240  //----------------------
241  Int_t StFcsPi0ReconstructionMaker::Make()
242  {
243 
244  StEvent* event = (StEvent*)GetInputDS("StEvent");
245  if(!event) {LOG_ERROR << "StFcsPi0ReconstructionMaker::Make did not find StEvent"<<endm; return kStErr;}
246  mFcsColl = event->fcsCollection();
247  if(!mFcsColl) {LOG_ERROR << "StFcsPi0ReconstructionMaker::Make did not find StEvent->StFcsCollection"<<endm; return kStErr;}
248 
249  h0->Fill(1); //count total number of event
250 
251  if(mNAccepted < mMaxEvents){
252  if(mFilter==1 && mFcsColl->numberOfHits(0)+ mFcsColl->numberOfHits(1)+ mFcsColl->numberOfHits(2)+ mFcsColl->numberOfHits(3)==0) return kStOK;
253  mNAccepted++;
254  int n_EcalMult=0;
255  int n_EcalClustMult=0;
256  int n_EcalPointMult=0;
257  int n_EcalClust_cut=0;
258  int n_EcalPoint_cut=0;
259  int n_Ecal_cut=0;
260 
261  for(int det=0; det<kFcsNDet; det++)
262  {
263  int check_ecal= mFcsDb->ecalHcalPres(det);
264  if (check_ecal!=0) continue;
265 
266  // get Ecal Mult, Ecal Cluster Mult
267  StSPtrVecFcsHit& hits = mFcsColl->hits(det);
268  int nh=mFcsColl->numberOfHits(det);
269  StSPtrVecFcsCluster& clusters = mFcsColl->clusters(det);
270  int nc=mFcsColl->numberOfClusters(det);
271  StSPtrVecFcsPoint& points = mFcsColl->points(det);
272  int np=mFcsColl->numberOfPoints(det);
273 
274  //count multiplicity
275  if ((det==0) || (det==1))
276  {
277  n_EcalMult=n_EcalMult + nh;
278  n_EcalClustMult=n_EcalClustMult + nc;
279  n_EcalPointMult=n_EcalPointMult + np;
280  };
281 
282  for (int i=0;i<nh;i++)
283  {
284  StFcsHit* hit=hits[i];
285  unsigned short hit_id=hit->id();
286  float hit_energy=hit->energy();
287  h1list_Etower[hit_id]->Fill(hit_energy); //fill in energy spectrum for tower
288  if (hit_energy>1) n_Ecal_cut++;//count Ecal Multiplicity when hit energy >1
289  }
290  for (int i=0;i<nc;i++)
291  {
292  StFcsCluster* clu=clusters[i];
293  float clu_energy=clu->energy();
294  if (clu_energy>1) n_EcalClust_cut++;//count Ecal cluster Multiplicity when cluster energy >1
295 
296  }
297  for (int i=0;i<np;i++)
298  {
299  StFcsPoint* pnt=points[i];
300  float pnt_energy=pnt->energy();
301  if (pnt_energy>1) n_EcalPoint_cut++;
302  }
303  }
304 
305  //TOF Mult
306  StTriggerData* trg=0;
307  StFcsRawDaqReader* fcsraw=(StFcsRawDaqReader*)GetMaker("daqReader");
308  StEvent* event= (StEvent*)GetInputDS("StEvent");
309  if(fcsraw)
310  {
311  trg = fcsraw->trgdata();
312  if(!trg){
313  LOG_INFO << "Canot find Trigger Data from StFcsRawDaqReader" << endm;
314  }
315  } else if(event){
316  trg=event->StEvent::triggerData();
317  if(!trg){
318  LOG_INFO << "Canot find Trigger Data from StEvent" << endm;
319  }
320  }
321 
322 
323  int tofmult = trg->tofMultiplicity();
324 
325  h2_EcalMult_vs_TOFMult->Fill(tofmult,n_EcalMult);
326  h2_EcalMult_vs_TOFMult_E_1->Fill(tofmult,n_Ecal_cut);
327 
328  if ((tofmult>TofMult_cut)||(tofmult<=2)) return kStOK;//Tof Mult cut
329 
330  if (n_Ecal_cut>EcalMult_cut) return kStOk;//Ecal Mult cut
331 
332  float h1_inv_mass_cluster_best=-1;
333  float h1_Zgg_cluster_best=-1;
334  float h1_opening_angle_cluster_best=-1;
335  float h1_two_cluster_energy_aftercut_best=-1;
336  float h1_dgg_cluster_best=-1;
337  float e1_cluster_best=-1;
338  float e2_cluster_best=-1;
339  float energy_largest_cluster=-1;
340  float clux1_best=-1;
341  float cluy1_best=-1;
342  float clux2_best=-1;
343  float cluy2_best=-1;
344  unsigned short best_tower_id1_cluster=0;
345  unsigned short best_tower_id2_cluster=0;
346 
347  float h1_inv_mass_point_best=-1;
348  float h1_Zgg_point_best=-1;
349  float h1_opening_angle_point_best=-1;
350  float h1_two_point_energy_aftercut_best=-1;
351  float h1_dgg_point_best=-1;
352  float e1_point_best=-1;
353  float e2_point_best=-1;
354  float energy_largest_point=-1;
355  float pntx1_best=-1;
356  float pnty1_best=-1;
357  float pntx2_best=-1;
358  float pnty2_best=-1;
359  float pntx1_rPosition_best=-1;
360  float pnty1_rPosition_best=-1;
361  float pntx2_rPosition_best=-1;
362  float pnty2_rPosition_best=-1;
363  unsigned short best_tower_id1_point=0;
364  unsigned short best_tower_id2_point=0;
365 
366  for(int det=0; det<kFcsNDet; det++)
367  {
368  int check_ecal= mFcsDb->ecalHcalPres(det);
369  if (check_ecal!=0) continue;
370  StSPtrVecFcsCluster& clusters = mFcsColl->clusters(det);
371  int nc=mFcsColl->numberOfClusters(det);
372  StSPtrVecFcsPoint& points = mFcsColl->points(det);
373  int np=mFcsColl->numberOfPoints(det);
374 
375  //no cut (cluster finder)
376 
377  for (int i=0;i<nc;i++)
378  {
379  StFcsCluster* clu=clusters[i];
380  float clu_energy=clu->energy();
381  h1_each_cluster_energy_nocut->Fill(clu_energy);
382 
383  if (i==nc-1) continue;
384  for (int j=i+1;j<nc;j++)
385  {
386  StFcsCluster* cluj=clusters[j];
387  float cluj_energy=cluj->energy();
388  h1_two_cluster_energy_nocut->Fill(clu_energy+cluj_energy);
389  }
390  }
391 
392  //no cut (point maker)
393  for (int i=0;i<np;i++)
394  {
395  StFcsPoint* pnt=points[i];
396  float pnt_energy=pnt->energy();
397  h1_each_point_energy_nocut->Fill(pnt_energy);
398 
399  if (i==np-1) continue;
400  for (int j=i+1;j<np;j++)
401  {
402  StFcsPoint* pntj=points[j];
403  float pntj_energy=pntj->energy();
404  h1_two_point_energy_nocut->Fill(pnt_energy+pntj_energy);
405  }
406  }
407 
408  //start cut (cluster finder)
409  for (int i=0; i<(nc-1); i++)
410  {
411 
412  StFcsCluster* clu=clusters[i];
413  float clu_x=clu->x();
414  float clu_y=clu->y();
415  float clu_energy=clu->energy();
416  if (clu_energy<1) continue; //cluster energy cut
417  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(det,clu_x,clu_y);
418  StLorentzVectorD p = mFcsDb->getLorentzVector((xyz),clu_energy,0);
419 
420 
421  for (int j=i+1;j<nc;j++)
422  {
423 
424  StFcsCluster* cluj=clusters[j];
425  float cluj_energy=cluj->energy();
426 
427  if ((cluj->energy())<1) continue;
428  float zgg=(fabs(clu_energy-cluj_energy))/(clu_energy+cluj_energy);
429  if (zgg>=0.7) continue;//Zgg cut
430  if (((clu_energy+cluj_energy)<8)||((clu_energy+cluj_energy)>20)) continue;//E1 + E2 cut
431 
432  StThreeVectorD xyzj = mFcsDb->getStarXYZfromColumnRow(det,cluj->x(),cluj->y());
433  StLorentzVectorD pj = mFcsDb->getLorentzVector((xyzj),cluj->energy(),0);
434 
435  //choose the best pair of cluster
436  if ((clu_energy+cluj_energy)>energy_largest_cluster)
437  {
438  energy_largest_cluster=clu_energy+cluj_energy;
439  h1_Zgg_cluster_best=zgg;
440  h1_inv_mass_cluster_best=(p+pj).m();
441  h1_two_cluster_energy_aftercut_best=clu_energy+cluj_energy;
442  h1_opening_angle_cluster_best=(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]))));
443  h1_dgg_cluster_best=(sqrt((xyz[0]-xyzj[0])*(xyz[0]-xyzj[0]) + (xyz[1]-xyzj[1])*(xyz[1]-xyzj[1])));
444  e1_cluster_best=clu_energy;
445  e2_cluster_best=cluj_energy;
446  clux1_best=clu->x();
447  cluy1_best=clu->y();
448  clux2_best=cluj->x();
449  cluy2_best=cluj->y();
450 
451  StPtrVecFcsHit& clui_hits= clu->hits();
452  StPtrVecFcsHit& cluj_hits= cluj->hits();
453  int clui_nTowers= clu->nTowers();
454  int cluj_nTowers= cluj->nTowers();
455  float max_energy_tower = -1;
456  unsigned short tower_id = 0;
457  for (int k=0;k<clui_nTowers;k++)
458  {
459  if (clui_hits[k]->energy()>max_energy_tower)
460  { max_energy_tower=clui_hits[k]->energy();
461  tower_id = clui_hits[k]->id(); }
462  }
463  best_tower_id1_cluster=tower_id;
464  for (int k=0;k<cluj_nTowers;k++)
465  {
466  if ((cluj_hits[k]->energy())>max_energy_tower)
467  { max_energy_tower=cluj_hits[k]->energy();
468  tower_id = cluj_hits[k]->id(); }
469  }
470  best_tower_id2_cluster=tower_id;
471 
472  }
473  }
474  }
475 
476 
477  //start cut(point)
478  for (int i=0; i<(np-1); i++)
479  {
480 
481  StFcsPoint* pnt=points[i];
482  float pnt_x=pnt->x();
483  float pnt_y=pnt->y();
484  float pnt_energy=pnt->energy();
485  if (pnt_energy<1) continue; //point energy cut
486  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(det,pnt_x,pnt_y);
487  StLorentzVectorD p = mFcsDb->getLorentzVector((xyz),pnt_energy,0);
488 
489  for (int j=i+1;j<np;j++)
490  {
491 
492  StFcsPoint* pntj=points[j];
493  float pntj_energy=pntj->energy();
494 
495  if ((pntj->energy())<1) continue;
496  float zgg=(fabs(pnt_energy-pntj_energy))/(pnt_energy+pntj_energy);
497  if (zgg>=0.7) continue;
498  if (((pnt_energy+pntj_energy)<8)||((pnt_energy+pntj_energy)>20)) continue;
499 
500  StThreeVectorD xyzj = mFcsDb->getStarXYZfromColumnRow(det,pntj->x(),pntj->y());
501  StLorentzVectorD pj = mFcsDb->getLorentzVector((xyzj),pntj->energy(),0);
502 
503  //choose the best pair of point
504  if ((pnt_energy+pntj_energy)>energy_largest_point)
505  {
506  energy_largest_point=pnt_energy+pntj_energy;
507  h1_Zgg_point_best=zgg;
508  h1_inv_mass_point_best=(p+pj).m();
509  h1_two_point_energy_aftercut_best=pnt_energy+pntj_energy;
510  h1_opening_angle_point_best=(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]))));
511  h1_dgg_point_best=(sqrt((xyz[0]-xyzj[0])*(xyz[0]-xyzj[0]) + (xyz[1]-xyzj[1])*(xyz[1]-xyzj[1])));
512  e1_point_best=pnt_energy;
513  e2_point_best=pntj_energy;
514  pntx1_best=pnt->x();
515  pnty1_best=pnt->y();
516  pntx2_best=pntj->x();
517  pnty2_best=pntj->y();
518  pntx1_rPosition_best=pntx1_best-int(pntx1_best);
519  pnty1_rPosition_best=pnty1_best-int(pnty1_best);
520  pntx2_rPosition_best=pntx2_best-int(pntx2_best);
521  pnty2_rPosition_best=pnty2_best-int(pnty2_best);
522 
523  StFcsCluster* pnt_clu=pnt->cluster();
524  StFcsCluster* pntj_clu=pntj->cluster();
525 
526  StPtrVecFcsHit& pnti_hits= pnt_clu->hits(); StPtrVecFcsHit& pntj_hits= pntj_clu->hits();
527  int pnti_nTowers= pnt_clu->nTowers();
528  int pntj_nTowers= pntj_clu->nTowers();
529  float max_energy_tower = -1;
530  unsigned short tower_id = 0;
531  for (int k=0;k<pnti_nTowers;k++)
532  {
533  if (pnti_hits[k]->energy()>max_energy_tower)
534  { max_energy_tower=pnti_hits[k]->energy();
535  tower_id = pnti_hits[k]->id(); }
536  }
537  best_tower_id1_point=tower_id;
538  for (int k=0;k<pntj_nTowers;k++)
539  {
540  if ((pntj_hits[k]->energy())>max_energy_tower)
541  { max_energy_tower=pntj_hits[k]->energy();
542  tower_id = pntj_hits[k]->id(); }
543  }
544  best_tower_id2_point=tower_id;
545 
546  }
547  }
548 
549 
550  }
551 
552  }
553 
554  //output
555  if (h1_inv_mass_cluster_best>-1)
556  {
557 
558  h1_inv_mass_cluster->Fill(h1_inv_mass_cluster_best);
559  h1_Zgg_cluster->Fill(h1_Zgg_cluster_best);
560  h1_opening_angle_cluster->Fill(h1_opening_angle_cluster_best);
561  h1_each_cluster_energy_aftercut->Fill(e1_cluster_best);
562  h1_each_cluster_energy_aftercut->Fill(e2_cluster_best);
563  h1_two_cluster_energy_aftercut->Fill(h1_two_cluster_energy_aftercut_best);
564  h1_dgg_cluster->Fill(h1_dgg_cluster_best);
565  h2_EcalClustMult_vs_TOFMult->Fill(tofmult,n_EcalClustMult);
566  h2_EcalClustMult_vs_TOFMult_E_1->Fill(tofmult,n_EcalClust_cut);
567  h2_dgg_cluster_vs_E1pE2_cluster->Fill(h1_two_cluster_energy_aftercut_best,h1_dgg_cluster_best);
568  h2_mass_cluster_vs_Zgg_cluster->Fill(h1_Zgg_cluster_best,h1_inv_mass_cluster_best);
569  h2_cluster_position->Fill(clux1_best,cluy1_best);
570  h2_cluster_position->Fill(clux2_best,cluy2_best);
571  h2_mass_cluster_vs_dgg_cluster->Fill(h1_dgg_cluster_best,h1_inv_mass_cluster_best);
572  h1list_mass_cluster_by_tower[best_tower_id1_cluster]->Fill(h1_inv_mass_cluster_best);
573  h1list_mass_cluster_by_tower[best_tower_id2_cluster]->Fill(h1_inv_mass_cluster_best);
574 
575  }
576 
577  //single event output for point
578  if (h1_inv_mass_point_best>-1)
579  {
580 
581  h1_inv_mass_point->Fill(h1_inv_mass_point_best);
582  h1_Zgg_point->Fill(h1_Zgg_point_best);
583  h1_opening_angle_point->Fill(h1_opening_angle_point_best);
584  h1_each_point_energy_aftercut->Fill(e1_point_best);
585  h1_each_point_energy_aftercut->Fill(e2_point_best);
586  h1_two_point_energy_aftercut->Fill(h1_two_point_energy_aftercut_best);
587  h1_dgg_point->Fill(h1_dgg_point_best);
588  h1_x_rPosition_point->Fill(pntx1_rPosition_best);
589  h1_y_rPosition_point->Fill(pnty1_rPosition_best);
590  h1_x_rPosition_point->Fill(pntx2_rPosition_best);
591  h1_y_rPosition_point->Fill(pnty2_rPosition_best);
592 
593  h2_EcalPointMult_vs_TOFMult->Fill(tofmult,n_EcalPointMult);
594  h2_EcalPointMult_vs_TOFMult_E_1->Fill(tofmult,n_EcalPoint_cut);
595  h2_dgg_point_vs_E1pE2_point->Fill(h1_two_point_energy_aftercut_best,h1_dgg_point_best);
596  h2_mass_point_vs_Zgg_point->Fill(h1_Zgg_point_best,h1_inv_mass_point_best);
597  h2_point_position->Fill(pntx1_best,pnty1_best);
598  h2_point_position->Fill(pntx2_best,pnty2_best);
599  h2_mass_point_vs_dgg_point->Fill(h1_dgg_point_best,h1_inv_mass_point_best);
600 
601  h1list_mass_point_by_tower[best_tower_id1_point]->Fill(h1_inv_mass_point_best);
602  h1list_mass_point_by_tower[best_tower_id2_point]->Fill(h1_inv_mass_point_best);
603 
604  }
605  }
606 
607  return kStOK;
608 }
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:853
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:821
StFcsDb::ecalHcalPres
int ecalHcalPres(int det) const
Ecal North=0, Ecal South=1, Hcal North=2, Hcal South=3, Pres=4/5.
Definition: StFcsDb.cxx:381
kStOK
Definition: Stypes.h:40
kStErr
Definition: Stypes.h:44
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5