Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JEventProcessor_pi0fcalskim.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: JEventProcessor_pi0fcalskim.cc
4 // Created: Mon Dec 1 14:57:11 EST 2014
5 // Creator: shepherd (on Linux ifarm1101 2.6.32-220.7.1.el6.x86_64 x86_64)
6 //
7 
8 #include <math.h>
9 
11 using namespace jana;
12 
13 // Routine used to create our JEventProcessor
14 #include "JANA/JApplication.h"
15 #include <TLorentzVector.h>
16 #include "TMath.h"
17 #include "JANA/JApplication.h"
18 #include "DANA/DApplication.h"
19 #include "FCAL/DFCALShower.h"
20 #include "FCAL/DFCALCluster.h"
21 #include "FCAL/DFCALHit.h"
23 #include "PID/DVertex.h"
24 #include "PID/DEventRFBunch.h"
25 
26 #include "GlueX.h"
27 #include <vector>
28 #include <deque>
29 #include <string>
30 #include <iostream>
31 #include <algorithm>
32 #include <stdio.h>
33 #include <stdlib.h>
34 
35 extern "C"{
36  void InitPlugin(JApplication *app){
37  InitJANAPlugin(app);
38  app->AddProcessor(new JEventProcessor_pi0fcalskim());
39  }
40 } // "C"
41 
42 
43 //------------------
44 // JEventProcessor_pi0fcalskim (Constructor)
45 //------------------
47 {
48 
49  WRITE_EVIO = 1;
50 /*
51  MIN_MASS = 0.03; // GeV
52  MAX_MASS = 0.30; // GeV
53  MIN_E = 1.0; // GeV (photon energy cut)
54  MIN_R = 20; // cm (cluster distance to beam line)
55  MAX_DT = 10; // ns (cluster time diff. cut)
56  MAX_ETOT = 12; // GeV (max total FCAL energy)
57  MIN_BLOCKS = 2; // minumum blocks per cluster
58 
59  WRITE_ROOT = 0;
60  WRITE_EVIO = 1;
61 
62  gPARMS->SetDefaultParameter( "PI0FCALSKIM:MIN_MASS", MIN_MASS );
63  gPARMS->SetDefaultParameter( "PI0FCALSKIM:MAX_MASS", MAX_MASS );
64  gPARMS->SetDefaultParameter( "PI0FCALSKIM:MIN_E", MIN_E );
65  gPARMS->SetDefaultParameter( "PI0FCALSKIM:MIN_R", MIN_R );
66  gPARMS->SetDefaultParameter( "PI0FCALSKIM:MAX_DT", MAX_DT );
67  gPARMS->SetDefaultParameter( "PI0FCALSKIM:MAX_ETOT", MAX_ETOT );
68  gPARMS->SetDefaultParameter( "PI0FCALSKIM:MIN_BLOCKS", MIN_BLOCKS );
69  gPARMS->SetDefaultParameter( "PI0FCALSKIM:WRITE_ROOT", WRITE_ROOT );
70  gPARMS->SetDefaultParameter( "PI0FCALSKIM:WRITE_EVIO", WRITE_EVIO );
71  */
72 }
73 
74 //------------------
75 // ~JEventProcessor_pi0fcalskim (Destructor)
76 //------------------
78 {
79 
80 }
81 
82 //------------------
83 // init
84 //------------------
86 {
87  num_epics_events = 0;
88 /*
89  if( ! ( WRITE_ROOT || WRITE_EVIO ) ){
90 
91  cerr << "No output mechanism has been specified." << endl;
92  return UNRECOVERABLE_ERROR;
93  }
94 
95  if( WRITE_ROOT ){
96 
97  japp->RootWriteLock();
98 
99  m_tree = new TTree( "cluster", "Cluster Tree for Pi0 Calibration" );
100  m_tree->Branch( "nClus", &m_nClus, "nClus/I" );
101  m_tree->Branch( "hit0", m_hit0, "hit0[nClus]/I" );
102  m_tree->Branch( "px", m_px, "px[nClus]/F" );
103  m_tree->Branch( "py", m_py, "py[nClus]/F" );
104  m_tree->Branch( "pz", m_pz, "pz[nClus]/F" );
105 
106  m_tree->Branch( "nHit", &m_nHit, "nHit/I" );
107  m_tree->Branch( "chan", m_chan, "chan[nHit]/I" );
108  m_tree->Branch( "e", m_e, "e[nHit]/F" );
109 
110  japp->RootUnLock();
111  }
112 */
113  return NOERROR;
114 }
115 
116 //------------------
117 // brun
118 //------------------
119 jerror_t JEventProcessor_pi0fcalskim::brun(JEventLoop *eventLoop, int32_t runnumber)
120 {
121  return NOERROR;
122 }
123 
124 //------------------
125 // evnt
126 //------------------
127 jerror_t JEventProcessor_pi0fcalskim::evnt(JEventLoop *loop, uint64_t eventnumber)
128 {
129 
130  vector< const DFCALShower* > locFCALShowers;
131  vector< const DVertex* > kinfitVertex;
132  loop->Get(locFCALShowers);
133  loop->Get(kinfitVertex);
134 
135  vector< const DTrackTimeBased* > locTrackTimeBased;
136  loop->Get(locTrackTimeBased);
137 
138  vector < const DFCALShower * > matchedShowers;
139 
140  const DEventWriterEVIO* locEventWriterEVIO = NULL;
141  loop->GetSingle(locEventWriterEVIO);
142 
143  // always write out BOR events
144  if(loop->GetJEvent().GetStatusBit(kSTATUS_BOR_EVENT)) {
145  //jout << "Found BOR!" << endl;
146  locEventWriterEVIO->Write_EVIOEvent( loop, "pi0fcalskim" );
147  return NOERROR;
148  }
149 
150  // write out the first few EPICS events to save run number & other meta info
151  if(loop->GetJEvent().GetStatusBit(kSTATUS_EPICS_EVENT) && (num_epics_events<5)) {
152  //jout << "Found EPICS!" << endl;
153  locEventWriterEVIO->Write_EVIOEvent( loop, "pi0fcalskim" );
154  num_epics_events++;
155  return NOERROR;
156  }
157 
158  vector< const JObject* > locObjectsToSave;
159 
160  bool Candidate = false;
161 
162  Double_t kinfitVertexX = 0.0, kinfitVertexY = 0.0, kinfitVertexZ = 0.0;
163 
164  for (unsigned int i = 0 ; i < kinfitVertex.size(); i++)
165  {
166  if(i==0)
167  locObjectsToSave.push_back(static_cast<const JObject *>(kinfitVertex[0]));
168 
169  kinfitVertexX = kinfitVertex[i]->dSpacetimeVertex.X();
170  kinfitVertexY = kinfitVertex[i]->dSpacetimeVertex.Y();
171  kinfitVertexZ = kinfitVertex[i]->dSpacetimeVertex.Z();
172  }
173 
174  vector<const DEventRFBunch*> locEventRFBunches;
175  loop->Get(locEventRFBunches);
176  if(locEventRFBunches.size() > 0) {
177  locObjectsToSave.push_back(static_cast<const JObject *>(locEventRFBunches[0]));
178  }
179 
180  DVector3 norm(0.0,0.0,-1);
181  DVector3 pos,mom;
182  // Double_t radius = 0;
183  //japp->RootWriteLock();
184  //Double_t p;
185  for (unsigned int i=0; i < locTrackTimeBased.size() ; ++i){
186  vector<DTrackFitter::Extrapolation_t>extrapolations=locTrackTimeBased[i]->extrapolations.at(SYS_FCAL);
187  if (extrapolations.size()==0) continue;
188 
189  for (unsigned int j=0; j< locFCALShowers.size(); ++j){
190 
191  Double_t x = locFCALShowers[j]->getPosition().X();
192  Double_t y = locFCALShowers[j]->getPosition().Y();
193  // Double_t z = locFCALShowers[j]->getPosition().Z() ;
194  //cout << "Z: " << z << endl;
195  //DVector3 pos_FCAL(x,y,625.406);
196  //for LH2 target
197  //DVector3 pos_FCAL(0,0,625.406);
198 
199  DVector3 pos_FCAL(0,0,638);
200  //at the end of the start counter; use this fall for fall '15 data
201  // DVector3 pos_FCAL(0,0,692);
202  //DVector3 pos_FCAL(0.0,0.0,650);
203  //std::cout<<"i: "<< i<< " j: "<< j << " z: "<<z<< endl;
204  // if (locTrackTimeBased[i]->rt->GetIntersectionWithPlane(pos_FCAL,norm,pos,mom)==NOERROR)
205  pos=extrapolations[0].position;
206 
207  // Double_t dX= TMath::Abs(pos.X() - x);
208  // Double_t dY= TMath::Abs(pos.Y() - y);
209  // Double_t dZ= TMath::Abs(pos.Z() - z);
210  Double_t trkmass = locTrackTimeBased[i]->mass();
211  Double_t FOM = TMath::Prob(locTrackTimeBased[i]->chisq, locTrackTimeBased[i]->Ndof);
212  // radius = sqrt(pos.X()*pos.X() + pos.Y()*pos.Y());
213  // Double_t Eshwr = locFCALShowers[j]->getEnergy();
214  // p = locTrackTimeBased[i]->momentum().Mag();
215  // cout<<"p: "<<p<<endl;
216  // Double_t dZ = TMath::Abs(pos.Z() - z);
217  Double_t dRho = sqrt(((pos.X() - x)*(pos.X() - x)) + ((pos.Y() - y)* (pos.Y() - y)));
218  // std::cout<<"i: "<< i<< " j: "<< j << " dRho " <<dRho << endl;
219  //if(dX < 20 && dY < 20 && trkmass < 0.15 && dRho < 50 && FOM > 0.01) {
220  if(trkmass < 0.15 && dRho < 5 && FOM > 0.01 ) {
221  matchedShowers.push_back(locFCALShowers[j]);
222  // matchedTracks.push_back(locTrackTimeBased[i]);
223  // printf ("Matched event=%d, i=%d, j=%d, p=%f, Ztrk=%f Zshr=%f, Xtrk=%f, Xshr=%f, Ytrk=%f, Yshr=%f\n",locEventNumber,i,j,p,
224  // pos.Z(),z,pos.X(),x,pos.Y(),y);
225  // break;
226 
227  }
228 
229  }
230  }
231 
232  for(unsigned int i=0; i<locFCALShowers.size(); i++)
233  {
234  if (find(matchedShowers.begin(), matchedShowers.end(),locFCALShowers[i]) != matchedShowers.end()) continue;
235 
236  const DFCALShower *s1 = locFCALShowers[i];
237 
238  vector<const DFCALCluster*> associated_clusters1;
239 
240  s1->Get(associated_clusters1);
241  Double_t dx1 = s1->getPosition().X() - kinfitVertexX;
242  Double_t dy1 = s1->getPosition().Y() - kinfitVertexY;
243  Double_t dz1 = s1->getPosition().Z() - kinfitVertexZ;
244  Double_t R1 = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
245  Double_t E1 = s1->getEnergy();
246  Double_t t1 = s1->getTime();
247  TLorentzVector sh1_p(E1*dx1/R1, E1*dy1/R1, E1*dz1/R1, E1);
248 
249  for(unsigned int j=i+1; j<locFCALShowers.size(); j++)
250  {
251  const DFCALShower *s2 = locFCALShowers[j];
252  if (find(matchedShowers.begin(), matchedShowers.end(),s2) != matchedShowers.end()) continue;
253 
254  vector<const DFCALCluster*> associated_clusters2;
255  s2->Get(associated_clusters2);
256  Double_t dx2 = s2->getPosition().X() - kinfitVertexX;
257  Double_t dy2 = s2->getPosition().Y() - kinfitVertexY;
258  Double_t dz2 = s2->getPosition().Z() - kinfitVertexZ;
259  Double_t R2 = sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
260  Double_t E2 = s2->getEnergy();
261  Double_t t2 = s2->getTime();
262 
263  TLorentzVector sh2_p(E2*dx2/R2, E2*dy2/R2, E2*dz2/R2, E2);
264  TLorentzVector ptot = sh1_p+sh2_p;
265  Double_t inv_mass = ptot.M();
266 
267  //Candidate |= (E1 > 0.5 && E2 > 0.5 && s1->getPosition().Pt() > 20*k_cm && s2->getPosition().Pt() > 20*k_cm && (fabs (t1-t2) < 10) && (inv_mass<0.30) ) ;
268  Candidate |= (E1 > 0.5 && E2 > 0.5 && (fabs (t1-t2) < 10) && (inv_mass<0.30) ) ;
269 
270  //if(E1 > 0.5 && E2 > 0.5 && s1->getPosition().Pt() > 20*k_cm && s2->getPosition().Pt() > 20*k_cm && (fabs (t1-t2) < 10) && (inv_mass<0.30) ) {
271  if(E1 > 0.5 && E2 > 0.5 && (fabs (t1-t2) < 10) && (inv_mass<0.30) ) {
272  if(find(locObjectsToSave.begin(), locObjectsToSave.end(), locFCALShowers[i]) == locObjectsToSave.end())
273  locObjectsToSave.push_back(static_cast<const JObject *>(locFCALShowers[i]));
274  if(find(locObjectsToSave.begin(), locObjectsToSave.end(), locFCALShowers[j]) == locObjectsToSave.end())
275  locObjectsToSave.push_back(static_cast<const JObject *>(locFCALShowers[j]));
276  }
277  }
278  }
279 
280  if( Candidate ){
281 
282  if( WRITE_EVIO ){
283 
284  locEventWriterEVIO->Write_EVIOEvent( loop, "pi0fcalskim", locObjectsToSave );
285  }
286  }
287 
288 
289  /*
290  vector< const DFCALCluster* > clusterVec;
291  loop->Get( clusterVec );
292 
293  if( clusterVec.size() < 2 ) return NOERROR;
294 
295  bool hasCandidate = false;
296  double eTot = 0;
297 
298  for( vector< const DFCALCluster*>::const_iterator clus1Itr = clusterVec.begin();
299  clus1Itr != clusterVec.end(); ++clus1Itr ){
300 
301  eTot += (**clus1Itr).getEnergy();
302 
303  for( vector< const DFCALCluster*>::const_iterator clus2Itr = clus1Itr + 1;
304  clus2Itr != clusterVec.end(); ++clus2Itr ){
305 
306  const DFCALCluster& clusL =
307  ( (**clus1Itr).getEnergy() > (**clus2Itr).getEnergy() ?
308  (**clus2Itr) : (**clus1Itr) );
309 
310  const DFCALCluster& clusH =
311  ( (**clus1Itr).getEnergy() > (**clus2Itr).getEnergy() ?
312  (**clus1Itr) : (**clus2Itr) );
313 
314  double clusLX = clusL.getCentroid().X();
315  double clusLY = clusL.getCentroid().Y();
316  double rL = sqrt( clusLX * clusLX + clusLY * clusLY );
317  double eL = clusL.getEnergy();
318  double tL = clusL.getTime();
319  int nHitL = clusL.GetHits().size();
320 
321  double clusHX = clusH.getCentroid().X();
322  double clusHY = clusH.getCentroid().Y();
323  double rH = sqrt( clusHX * clusHX + clusHY * clusHY );
324  double eH = clusH.getEnergy();
325  double tH = clusH.getTime();
326  int nHitH = clusH.GetHits().size();
327 
328  DVector3 clusLMom = clusL.getCentroid();
329  clusLMom.SetMag( eL );
330 
331  DVector3 clusHMom = clusH.getCentroid();
332  clusHMom.SetMag( eH );
333 
334  double dt = fabs( tL - tH );
335 
336  DLorentzVector gamL( clusLMom, clusLMom.Mag() );
337  DLorentzVector gamH( clusHMom, clusHMom.Mag() );
338 
339  double mass = ( gamL + gamH ).M();
340 
341  √ hasCandidate |=
342  ( ( eL > MIN_E ) &&
343  ( dt < MAX_DT ) &&
344  ( rL > MIN_R ) && ( rH > MIN_R ) &&
345 √√ ( nHitL >= MIN_BLOCKS ) && ( nHitH >= MIN_BLOCKS ) &&
346  ( mass > MIN_MASS ) && ( mass < MAX_MASS ) );
347  }
348  }
349 
350  if( hasCandidate && ( eTot < MAX_ETOT ) ){
351 
352  if( WRITE_EVIO ){
353 
354  dEventWriterEVIO->Write_EVIOEvent( loop, "pi0fcalskim" );
355  }
356 
357  if( WRITE_ROOT ){
358 
359  japp->RootWriteLock();
360  writeClustersToRoot( clusterVec );
361  japp->RootUnLock();
362  }
363  }
364 */
365  return NOERROR;
366 }
367 
368 //------------------
369 // erun
370 //------------------
372 {
373  // This is called whenever the run number changes, before it is
374  // changed to give you a chance to clean up before processing
375  // events from the next run number.
376  return NOERROR;
377 }
378 
379 //------------------
380 // Fin
381 //------------------
383 {
384  // Called before program exit after event processing is finished.
385  return NOERROR;
386 }
387 
388 
389 /*
390 void
391 JEventProcessor_pi0fcalskim::writeClustersToRoot( const vector< const DFCALCluster* > clusVec ){
392 
393  // this code must run serially -- obtain a lock before
394  // entering this function
395 
396  m_nHit = 0;
397  m_nClus = 0;
398 
399  // hit and cluster indices
400  int& iH = m_nHit;
401  int& iC = m_nClus;
402 
403  for( vector< const DFCALCluster* >::const_iterator clusItr = clusVec.begin();
404  clusItr != clusVec.end(); ++clusItr ){
405 
406  // if we exceed max clusters abort writing for this event
407  if( iC == kMaxClus ) return;
408 
409  const DFCALCluster& clus = (**clusItr);
410 
411  if( ( clus.getCentroid().Perp() < 20*k_cm ) ||
412  ( clus.getEnergy() < 1*k_GeV ) ||
413  ( clus.GetHits().size() < 2 ) ) continue;
414 
415  DVector3 gamMom = clus.getCentroid();
416  gamMom.SetMag( clus.getEnergy() );
417 
418  m_hit0[iC] = iH;
419  m_px[iC] = gamMom.X();
420  m_py[iC] = gamMom.Y();
421  m_pz[iC] = gamMom.Z();
422 
423  const vector<DFCALCluster::DFCALClusterHit_t>& hits = clus.GetHits();
424 
425  for( unsigned int i = 0; i < hits.size(); ++i ){
426 
427  // if we exceed max hits abort this event and return
428  if( iH == kMaxHits ) return;
429 
430  m_chan[iH] = hits[i].ch;
431  m_e[iH] = hits[i].E;
432  ++iH;
433  }
434  ++iC;
435  }
436 
437  m_tree->Fill();
438 }
439 */
double getEnergy() const
Definition: DFCALShower.h:156
static double E1[100]
jerror_t init(void)
Called once at program start.
TVector3 DVector3
Definition: DVector3.h:14
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
double getTime() const
Definition: DFCALShower.h:161
#define y
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
jerror_t fini(void)
Called after last event of last event source has been processed.
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
InitPlugin_t InitPlugin
TLatex * t1
Definition: GlueX.h:22
double sqrt(double)
bool Write_EVIOEvent(JEventLoop *locEventLoop, string locOutputFileNameSubString) const
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
DVector3 getPosition() const
Definition: DFCALShower.h:151