Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DEventWriterREST.cc
Go to the documentation of this file.
1 #include "DEventWriterREST.h"
2 
3 #include <DANA/DApplication.h>
4 #include <JANA/JCalibration.h>
5 
7 {
8  // must be read/used entirely in "RESTWriter" lock
9  static int locNumEventWriterThreads = 0;
10  return locNumEventWriterThreads;
11 }
12 
13 map<string, pair<ofstream*, hddm_r::ostream*> >& DEventWriterREST::Get_RESTOutputFilePointers(void) const
14 {
15  // must be read/used entirely in "RESTWriter" lock
16  // cannot do individual file locks, because the map itself can be modified
17  static map<string, pair<ofstream*, hddm_r::ostream*> > locRESTOutputFilePointers;
18  return locRESTOutputFilePointers;
19 }
20 
21 DEventWriterREST::DEventWriterREST(JEventLoop* locEventLoop, string locOutputFileBaseName) : dOutputFileBaseName(locOutputFileBaseName)
22 {
23  japp->WriteLock("RESTWriter");
24  {
26  }
27  japp->Unlock("RESTWriter");
28 
29  HDDM_USE_COMPRESSION = true;
30  string locCompressionString = "Turn on/off compression of the output HDDM stream. Set to \"0\" to turn off (it's on by default)";
31  gPARMS->SetDefaultParameter("HDDM:USE_COMPRESSION", HDDM_USE_COMPRESSION, locCompressionString);
32 
34  string locIntegrityString = "Turn on/off automatic integrity checking on the output HDDM stream. Set to \"0\" to turn off (it's on by default)";
35  gPARMS->SetDefaultParameter("HDDM:USE_INTEGRITY_CHECKS", HDDM_USE_INTEGRITY_CHECKS, locIntegrityString);
36 
38  if(gPARMS->Exists("REST:DATAVERSIONSTRING"))
39  gPARMS->GetParameter("REST:DATAVERSIONSTRING", HDDM_DATA_VERSION_STRING);
40  else
41  gPARMS->SetDefaultParameter("REST:DATAVERSIONSTRING", HDDM_DATA_VERSION_STRING, "");
42 
43  REST_WRITE_DIRC_HITS = true;
44  gPARMS->SetDefaultParameter("REST:WRITE_DIRC_HITS", REST_WRITE_DIRC_HITS);
45 
47  // if we can get the calibration context from the DANA interface, then save this as well
48  DApplication *dapp = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
49  if (dapp) {
50  JEvent& event = locEventLoop->GetJEvent();
51  JCalibration *jcalib = dapp->GetJCalibration(event.GetRunNumber());
52  if (jcalib) {
53  CCDB_CONTEXT_STRING = jcalib->GetContext();
54  }
55  }
56 
57 }
58 
59 bool DEventWriterREST::Write_RESTEvent(JEventLoop* locEventLoop, string locOutputFileNameSubString) const
60 {
61  std::vector<const DMCReaction*> reactions;
62  locEventLoop->Get(reactions);
63 
64  std::vector<const DRFTime*> rftimes;
65  locEventLoop->Get(rftimes);
66 
67  std::vector<const DBeamPhoton*> locBeamPhotons;
68  locEventLoop->Get(locBeamPhotons);
69 
70  std::vector<const DBeamPhoton*> locBeamPhotons_TAGGEDMCGEN;
71  locEventLoop->Get(locBeamPhotons_TAGGEDMCGEN, "TAGGEDMCGEN");
72 
73  std::vector<const DFCALShower*> fcalshowers;
74  locEventLoop->Get(fcalshowers);
75 
76  std::vector<const DBCALShower*> bcalshowers;
77  locEventLoop->Get(bcalshowers);
78 
79  std::vector<const DTOFPoint*> tofpoints;
80  locEventLoop->Get(tofpoints);
81 
82  std::vector<const DSCHit*> starthits;
83  locEventLoop->Get(starthits);
84 
85  std::vector<const DTrackTimeBased*> tracks;
86  locEventLoop->Get(tracks);
87 
88  std::vector<const DDetectorMatches*> locDetectorMatches;
89  locEventLoop->Get(locDetectorMatches);
90 
91  std::vector<const DDIRCPmtHit*> locDIRCPmtHits;
92  locEventLoop->Get(locDIRCPmtHits);
93 
94  std::vector<const DTrigger*> locTriggers;
95  locEventLoop->Get(locTriggers);
96 
97  //Check to see if there are any objects to write out. If so, don't write out an empty event
98  bool locOutputDataPresentFlag = false;
99  if((!reactions.empty()) || (!locBeamPhotons.empty()) || (!tracks.empty()))
100  locOutputDataPresentFlag = true;
101  else if((!fcalshowers.empty()) || (!bcalshowers.empty()) || (!tofpoints.empty()) || (!starthits.empty()))
102  locOutputDataPresentFlag = true;
103  //don't need to check detector matches: no matches if none of the above objects
104  if(!locOutputDataPresentFlag)
105  return true; //had correct response to data
106 
107  string locOutputFileName = Get_OutputFileName(locOutputFileNameSubString);
108 
109  hddm_r::HDDM locRecord;
110  hddm_r::ReconstructedPhysicsEventList res = locRecord.addReconstructedPhysicsEvents(1);
111 
112  // load the run and event numbers
113  JEvent& event = locEventLoop->GetJEvent();
114  res().setRunNo(event.GetRunNumber());
115  //The REST type for this is int64_t, whereas the event type is uint64_t
116  //This copy is lazy: the last bit is lost. However, we should never need the last bit.
117  res().setEventNo(event.GetEventNumber());
118 
119  // push any DMCReaction objects to the output record
120  for (size_t i=0; i < reactions.size(); i++)
121  {
122  hddm_r::ReactionList rea = res().addReactions(1);
123  rea().setType(reactions[i]->type);
124  rea().setWeight(reactions[i]->weight);
125  rea().setEbeam(reactions[i]->beam.energy());
126  rea().setTargetType(reactions[i]->target.PID());
127 
128  if(i != 0)
129  break;
130 
131  std::vector<const DMCThrown*> throwns;
132  locEventLoop->Get(throwns);
133  hddm_r::VertexList ver = rea().getVertices();
134  DLorentzVector locPreviousX4(-9.9E9, -9.9E9, -9.9E9, -9.9E9);
135  for(size_t it=0; it < throwns.size(); ++it)
136  {
137  DLorentzVector locThrownX4(throwns[it]->position(), throwns[it]->time());
138  if((locThrownX4.T() != locPreviousX4.T()) || (locThrownX4.Vect() != locPreviousX4.Vect()))
139  {
140  //new vertex
141  ver = rea().addVertices(1);
142  hddm_r::OriginList ori = ver().addOrigins(1);
143  ori().setT(locThrownX4.T());
144  ori().setVx(locThrownX4.X());
145  ori().setVy(locThrownX4.Y());
146  ori().setVz(locThrownX4.Z());
147  locPreviousX4 = locThrownX4;
148  }
149 
150  hddm_r::ProductList pro = ver().addProducts(1);
151  pro().setId(throwns[it]->myid);
152  pro().setParentId(throwns[it]->parentid);
153  int pdgtype = throwns[it]->pdgtype;
154  if (pdgtype == 0)
155  pdgtype = PDGtype((Particle_t)throwns[it]->type);
156  pro().setPdgtype(pdgtype);
157  hddm_r::MomentumList mom = pro().addMomenta(1);
158  mom().setE(throwns[it]->energy());
159  mom().setPx(throwns[it]->px());
160  mom().setPy(throwns[it]->py());
161  mom().setPz(throwns[it]->pz());
162  }
163  }
164 
165  // push any DRFTime objects to the output record
166  for (size_t i=0; i < rftimes.size(); i++)
167  {
168  hddm_r::RFtimeList rf = res().addRFtimes(1);
169  rf().setTsync(rftimes[i]->dTime);
170  }
171 
172  // push any DBeamPhoton objects to the output record
173  for(size_t loc_i = 0; loc_i < locBeamPhotons.size(); ++loc_i)
174  {
175  if(locBeamPhotons[loc_i]->dSystem == SYS_TAGM)
176  {
177  hddm_r::TagmBeamPhotonList locTagmBeamPhotonList = res().addTagmBeamPhotons(1);
178  locTagmBeamPhotonList().setT(locBeamPhotons[loc_i]->time());
179  locTagmBeamPhotonList().setE(locBeamPhotons[loc_i]->energy());
180  }
181  else if(locBeamPhotons[loc_i]->dSystem == SYS_TAGH)
182  {
183  hddm_r::TaghBeamPhotonList locTaghBeamPhotonList = res().addTaghBeamPhotons(1);
184  locTaghBeamPhotonList().setT(locBeamPhotons[loc_i]->time());
185  locTaghBeamPhotonList().setE(locBeamPhotons[loc_i]->energy());
186  }
187  }
188  for(size_t loc_i = 0; loc_i < locBeamPhotons_TAGGEDMCGEN.size(); ++loc_i)
189  {
190  if(locBeamPhotons_TAGGEDMCGEN[loc_i]->dSystem == SYS_TAGM)
191  {
192  hddm_r::TagmBeamPhotonList locTagmBeamPhotonList = res().addTagmBeamPhotons(1);
193  locTagmBeamPhotonList().setJtag("TAGGEDMCGEN");
194  locTagmBeamPhotonList().setT(locBeamPhotons_TAGGEDMCGEN[loc_i]->time());
195  locTagmBeamPhotonList().setE(locBeamPhotons_TAGGEDMCGEN[loc_i]->energy());
196  }
197  else if(locBeamPhotons_TAGGEDMCGEN[loc_i]->dSystem == SYS_TAGH)
198  {
199  hddm_r::TaghBeamPhotonList locTaghBeamPhotonList = res().addTaghBeamPhotons(1);
200  locTaghBeamPhotonList().setJtag("TAGGEDMCGEN");
201  locTaghBeamPhotonList().setT(locBeamPhotons_TAGGEDMCGEN[loc_i]->time());
202  locTaghBeamPhotonList().setE(locBeamPhotons_TAGGEDMCGEN[loc_i]->energy());
203  }
204  }
205 
206  // push any DFCALShower objects to the output record
207  for (size_t i=0; i < fcalshowers.size(); i++)
208  {
209  hddm_r::FcalShowerList fcal = res().addFcalShowers(1);
210  DVector3 pos = fcalshowers[i]->getPosition();
211  fcal().setX(pos(0));
212  fcal().setY(pos(1));
213  fcal().setZ(pos(2));
214  fcal().setT(fcalshowers[i]->getTime());
215  fcal().setE(fcalshowers[i]->getEnergy());
216  fcal().setXerr(fcalshowers[i]->xErr());
217  fcal().setYerr(fcalshowers[i]->yErr());
218  fcal().setZerr(fcalshowers[i]->zErr());
219  fcal().setTerr(fcalshowers[i]->tErr());
220  fcal().setEerr(fcalshowers[i]->EErr());
221  fcal().setXycorr(fcalshowers[i]->XYcorr());
222  fcal().setXzcorr(fcalshowers[i]->XZcorr());
223  fcal().setYzcorr(fcalshowers[i]->YZcorr());
224  fcal().setEzcorr(fcalshowers[i]->EZcorr());
225  fcal().setTzcorr(fcalshowers[i]->ZTcorr());
226 
227  // further correlations (an extension of REST format so code is different.)
228  hddm_r::FcalCorrelationsList locFcalCorrelationsList = fcal().addFcalCorrelationses(1);
229  locFcalCorrelationsList().setEtcorr(fcalshowers[i]->ETcorr());
230  locFcalCorrelationsList().setExcorr(fcalshowers[i]->EXcorr());
231  locFcalCorrelationsList().setEycorr(fcalshowers[i]->EYcorr());
232  locFcalCorrelationsList().setTxcorr(fcalshowers[i]->XTcorr());
233  locFcalCorrelationsList().setTycorr(fcalshowers[i]->YTcorr());
234 
235  // add in classification based on MVA
236  //hddm_r::FcalShowerClassificationList locFcalShowerClassificationList = fcal().addFcalShowerClassifications(1);
237  // locFcalShowerClassificationList().setClassifierOuput(fcalshowers[i]->getClassifierOutput());
238 
239  // add in shower properties used for MVA algorithm, etc.
240 
241  hddm_r::FcalShowerPropertiesList locFcalShowerPropertiesList = fcal().addFcalShowerPropertiesList(1);
242  locFcalShowerPropertiesList().setDocaTrack(fcalshowers[i]->getDocaTrack());
243  locFcalShowerPropertiesList().setTimeTrack(fcalshowers[i]->getTimeTrack());
244  locFcalShowerPropertiesList().setSumU(fcalshowers[i]->getSumU());
245  locFcalShowerPropertiesList().setSumV(fcalshowers[i]->getSumV());
246  locFcalShowerPropertiesList().setE1E9(fcalshowers[i]->getE1E9());
247  locFcalShowerPropertiesList().setE9E25(fcalshowers[i]->getE9E25());
248  hddm_r::FcalShowerNBlocksList locFcalShowerNBlocksList = fcal().addFcalShowerNBlockses(1);
249  locFcalShowerNBlocksList().setNumBlocks(fcalshowers[i]->getNumBlocks());
250 
251  }
252 
253 
254  // push any DBCALShower objects to the output record
255  for (size_t i=0; i < bcalshowers.size(); i++)
256  {
257  hddm_r::BcalShowerList bcal = res().addBcalShowers(1);
258  DVector3 pos(bcalshowers[i]->x,bcalshowers[i]->y,bcalshowers[i]->z);
259  bcal().setX(bcalshowers[i]->x);
260  bcal().setY(bcalshowers[i]->y);
261  bcal().setZ(bcalshowers[i]->z);
262  bcal().setT(bcalshowers[i]->t);
263  bcal().setE(bcalshowers[i]->E);
264  bcal().setXerr(bcalshowers[i]->xErr());
265  bcal().setYerr(bcalshowers[i]->yErr());
266  bcal().setZerr(bcalshowers[i]->zErr());
267  bcal().setTerr(bcalshowers[i]->tErr());
268  bcal().setEerr(bcalshowers[i]->EErr());
269  bcal().setXycorr(bcalshowers[i]->XYcorr());
270  bcal().setXzcorr(bcalshowers[i]->XZcorr());
271  bcal().setYzcorr(bcalshowers[i]->YZcorr());
272  bcal().setEzcorr(bcalshowers[i]->EZcorr());
273  bcal().setTzcorr(bcalshowers[i]->ZTcorr());
274 
275  // further correlations (an extension of REST format so code is different.)
276  hddm_r::BcalCorrelationsList locBcalCorrelationsList = bcal().addBcalCorrelationses(1);
277  locBcalCorrelationsList().setEtcorr(bcalshowers[i]->ETcorr());
278  locBcalCorrelationsList().setExcorr(bcalshowers[i]->EXcorr());
279  locBcalCorrelationsList().setEycorr(bcalshowers[i]->EYcorr());
280  locBcalCorrelationsList().setTxcorr(bcalshowers[i]->XTcorr());
281  locBcalCorrelationsList().setTycorr(bcalshowers[i]->YTcorr());
282 
283  hddm_r::PreshowerList locPreShowerList = bcal().addPreshowers(1);
284  locPreShowerList().setPreshowerE(bcalshowers[i]->E_preshower);
285 
286  hddm_r::WidthList locWidthList = bcal().addWidths(1);
287  locWidthList().setSigLong(bcalshowers[i]->sigLong);
288  locWidthList().setSigTrans(bcalshowers[i]->sigTrans);
289  locWidthList().setSigTheta(bcalshowers[i]->sigTheta);
290 
291  //N_cell
292  hddm_r::BcalClusterList bcalcluster = bcal().addBcalClusters(1);
293  bcalcluster().setNcell(bcalshowers[i]->N_cell);
294 
295  hddm_r::BcalLayersList bcallayerdata = bcal().addBcalLayerses(1);
296  bcallayerdata().setE_L2(bcalshowers[i]->E_L2);
297  bcallayerdata().setE_L3(bcalshowers[i]->E_L3);
298  bcallayerdata().setE_L4(bcalshowers[i]->E_L4);
299  bcallayerdata().setRmsTime(bcalshowers[i]->rmsTime);
300  }
301 
302  // push any DTOFPoint objects to the output record
303  for (size_t i=0; i < tofpoints.size(); i++)
304  {
305  hddm_r::TofPointList tof = res().addTofPoints(1);
306  tof().setX(tofpoints[i]->pos(0));
307  tof().setY(tofpoints[i]->pos(1));
308  tof().setZ(tofpoints[i]->pos(2));
309  tof().setT(tofpoints[i]->t);
310  tof().setDE(tofpoints[i]->dE);
311 
312  //Status //Assume compiler optimizes multiplication
313  hddm_r::TofStatusList tofstatus = tof().addTofStatuses(1);
314  int locStatus = tofpoints[i]->dHorizontalBar + 45*tofpoints[i]->dVerticalBar;
315  locStatus += 45*45*tofpoints[i]->dHorizontalBarStatus + 45*45*4*tofpoints[i]->dVerticalBarStatus;
316  tofstatus().setStatus(locStatus);
317  }
318 
319  // push any DSCHit objects to the output record
320  for (size_t i=0; i < starthits.size(); i++)
321  {
322  hddm_r::StartHitList hit = res().addStartHits(1);
323  hit().setSector(starthits[i]->sector);
324  hit().setT(starthits[i]->t);
325  hit().setDE(starthits[i]->dE);
326  }
327 
329  // push any DDIRCPmtHit objects to the output record
330  for (size_t i=0; i < locDIRCPmtHits.size(); i++)
331  {
332  hddm_r::DircHitList hit = res().addDircHits(1);
333  hit().setCh(locDIRCPmtHits[i]->ch);
334  hit().setT(locDIRCPmtHits[i]->t);
335  hit().setTot(locDIRCPmtHits[i]->tot);
336  }
337  }
338 
339  // push any DTrackTimeBased objects to the output record
340  for (size_t i=0; i < tracks.size(); ++i)
341  {
342  hddm_r::ChargedTrackList tra = res().addChargedTracks(1);
343  tra().setCandidateId(tracks[i]->candidateid);
344  tra().setPtype(tracks[i]->PID());
345 
346  hddm_r::TrackFitList fit = tra().addTrackFits(1);
347  fit().setNdof(tracks[i]->Ndof);
348  fit().setChisq(tracks[i]->chisq);
349  fit().setX0(tracks[i]->x());
350  fit().setY0(tracks[i]->y());
351  fit().setZ0(tracks[i]->z());
352  fit().setPx(tracks[i]->px());
353  fit().setPy(tracks[i]->py());
354  fit().setPz(tracks[i]->pz());
355  fit().setT0(tracks[i]->time());
356  fit().setT0err(0.0);
357  fit().setT0det(SYS_CDC);
358 
359  const TMatrixFSym& errors = *(tracks[i]->TrackingErrorMatrix().get());
360  fit().setE11(errors(0,0));
361  fit().setE12(errors(0,1));
362  fit().setE13(errors(0,2));
363  fit().setE14(errors(0,3));
364  fit().setE15(errors(0,4));
365  fit().setE22(errors(1,1));
366  fit().setE23(errors(1,2));
367  fit().setE24(errors(1,3));
368  fit().setE25(errors(1,4));
369  fit().setE33(errors(2,2));
370  fit().setE34(errors(2,3));
371  fit().setE35(errors(2,4));
372  fit().setE44(errors(3,3));
373  fit().setE45(errors(3,4));
374  fit().setE55(errors(4,4));
375 
376  hddm_r::TrackFlagsList myflags = tra().addTrackFlagses(1);
377  myflags().setFlags(tracks[i]->flags);
378 
379  hddm_r::HitlayersList locHitLayers = tra().addHitlayerses(1);
380  locHitLayers().setCDCrings(tracks[i]->dCDCRings);
381  locHitLayers().setFDCplanes(tracks[i]->dFDCPlanes);
382 
383  vector<const DCDCTrackHit*> locCDCHits;
384  tracks[i]->Get(locCDCHits);
385  vector<const DFDCPseudo*> locFDCHits;
386  tracks[i]->Get(locFDCHits);
387 
388  hddm_r::ExpectedhitsList locExpectedHits = tra().addExpectedhitses(1);
389  //locExpectedHits().setMeasuredCDChits(locCDCHits.size());
390  //locExpectedHits().setMeasuredFDChits(locFDCHits.size());
391  locExpectedHits().setMeasuredCDChits(tracks[i]->measured_cdc_hits_on_track);
392  locExpectedHits().setMeasuredFDChits(tracks[i]->measured_fdc_hits_on_track);
393  //locExpectedHits().setMeasuredCDChits(tracks[i]->cdc_hit_usage.total_hits);
394  //locExpectedHits().setMeasuredFDChits(tracks[i]->fdc_hit_usage.total_hits);
395  locExpectedHits().setExpectedCDChits(tracks[i]->potential_cdc_hits_on_track);
396  locExpectedHits().setExpectedFDChits(tracks[i]->potential_fdc_hits_on_track);
397 
398  hddm_r::McmatchList locMCMatches = tra().addMcmatchs(1);
399  locMCMatches().setIthrown(tracks[i]->dMCThrownMatchMyID);
400  locMCMatches().setNumhitsmatch(tracks[i]->dNumHitsMatchedToThrown);
401 
402  if (tracks[i]->dNumHitsUsedFordEdx_FDC + tracks[i]->dNumHitsUsedFordEdx_CDC > 0)
403  {
404  hddm_r::DEdxDCList elo = tra().addDEdxDCs(1);
405  elo().setNsampleFDC(tracks[i]->dNumHitsUsedFordEdx_FDC);
406  elo().setNsampleCDC(tracks[i]->dNumHitsUsedFordEdx_CDC);
407  elo().setDxFDC(tracks[i]->ddx_FDC);
408  elo().setDxCDC(tracks[i]->ddx_CDC);
409  elo().setDEdxFDC(tracks[i]->ddEdx_FDC);
410  elo().setDEdxCDC(tracks[i]->ddEdx_CDC);
411  hddm_r::CDCAmpdEdxList elo2 = elo().addCDCAmpdEdxs(1);
412  elo2().setDxCDCAmp(tracks[i]->ddx_CDC_amp);
413  elo2().setDEdxCDCAmp(tracks[i]->ddEdx_CDC_amp);
414 
415  }
416  }
417 
418  // push any DTrigger objects to the output record
419  for (size_t i=0; i < locTriggers.size(); ++i)
420  {
421  hddm_r::TriggerList trigger = res().addTriggers(1);
422  trigger().setL1_trig_bits(Convert_UnsignedIntToSigned(locTriggers[i]->Get_L1TriggerBits()));
423  trigger().setL1_fp_trig_bits(Convert_UnsignedIntToSigned(locTriggers[i]->Get_L1FrontPanelTriggerBits()));
424  }
425 
426  // push any DDetectorMatches objects to the output record
427  for(size_t loc_i = 0; loc_i < locDetectorMatches.size(); ++loc_i)
428  {
429  hddm_r::DetectorMatchesList matches = res().addDetectorMatcheses(1);
430  for(size_t loc_j = 0; loc_j < tracks.size(); ++loc_j)
431  {
432  vector<shared_ptr<const DBCALShowerMatchParams>> locBCALShowerMatchParamsVector;
433  locDetectorMatches[loc_i]->Get_BCALMatchParams(tracks[loc_j], locBCALShowerMatchParamsVector);
434  for(size_t loc_k = 0; loc_k < locBCALShowerMatchParamsVector.size(); ++loc_k)
435  {
436  hddm_r::BcalMatchParamsList bcalList = matches().addBcalMatchParamses(1);
437  bcalList().setTrack(loc_j);
438 
439  const DBCALShower* locBCALShower = locBCALShowerMatchParamsVector[loc_k]->dBCALShower;
440  size_t locBCALindex = 0;
441  for(; locBCALindex < bcalshowers.size(); ++locBCALindex)
442  {
443  if(bcalshowers[locBCALindex] == locBCALShower)
444  break;
445  }
446  bcalList().setShower(locBCALindex);
447 
448  bcalList().setDeltaphi(locBCALShowerMatchParamsVector[loc_k]->dDeltaPhiToShower);
449  bcalList().setDeltaz(locBCALShowerMatchParamsVector[loc_k]->dDeltaZToShower);
450  bcalList().setDx(locBCALShowerMatchParamsVector[loc_k]->dx);
451  bcalList().setPathlength(locBCALShowerMatchParamsVector[loc_k]->dPathLength);
452  bcalList().setTflight(locBCALShowerMatchParamsVector[loc_k]->dFlightTime);
453  bcalList().setTflightvar(locBCALShowerMatchParamsVector[loc_k]->dFlightTimeVariance);
454  }
455 
456  vector<shared_ptr<const DFCALShowerMatchParams>> locFCALShowerMatchParamsVector;
457  locDetectorMatches[loc_i]->Get_FCALMatchParams(tracks[loc_j], locFCALShowerMatchParamsVector);
458  for (size_t loc_k = 0; loc_k < locFCALShowerMatchParamsVector.size(); ++loc_k)
459  {
460  hddm_r::FcalMatchParamsList fcalList = matches().addFcalMatchParamses(1);
461  fcalList().setTrack(loc_j);
462 
463  const DFCALShower* locFCALShower = locFCALShowerMatchParamsVector[loc_k]->dFCALShower;
464  size_t locFCALindex = 0;
465  for(; locFCALindex < fcalshowers.size(); ++locFCALindex)
466  {
467  if(fcalshowers[locFCALindex] == locFCALShower)
468  break;
469  }
470  fcalList().setShower(locFCALindex);
471 
472  fcalList().setDoca(locFCALShowerMatchParamsVector[loc_k]->dDOCAToShower);
473  fcalList().setDx(locFCALShowerMatchParamsVector[loc_k]->dx);
474  fcalList().setPathlength(locFCALShowerMatchParamsVector[loc_k]->dPathLength);
475  fcalList().setTflight(locFCALShowerMatchParamsVector[loc_k]->dFlightTime);
476  fcalList().setTflightvar(locFCALShowerMatchParamsVector[loc_k]->dFlightTimeVariance);
477  }
478 
479  vector<shared_ptr<const DTOFHitMatchParams>> locTOFHitMatchParamsVector;
480  locDetectorMatches[loc_i]->Get_TOFMatchParams(tracks[loc_j], locTOFHitMatchParamsVector);
481  for(size_t loc_k = 0; loc_k < locTOFHitMatchParamsVector.size(); ++loc_k)
482  {
483  hddm_r::TofMatchParamsList tofList = matches().addTofMatchParamses(1);
484  tofList().setTrack(loc_j);
485 
486  size_t locTOFindex = 0;
487  for(; locTOFindex < tofpoints.size(); ++locTOFindex)
488  {
489  if(tofpoints[locTOFindex] == locTOFHitMatchParamsVector[loc_k]->dTOFPoint)
490  break;
491  }
492  tofList().setHit(locTOFindex);
493 
494  tofList().setThit(locTOFHitMatchParamsVector[loc_k]->dHitTime);
495  tofList().setThitvar(locTOFHitMatchParamsVector[loc_k]->dHitTimeVariance);
496  tofList().setEhit(locTOFHitMatchParamsVector[loc_k]->dHitEnergy);
497 
498  tofList().setDEdx(locTOFHitMatchParamsVector[loc_k]->dEdx);
499  tofList().setPathlength(locTOFHitMatchParamsVector[loc_k]->dPathLength);
500  tofList().setTflight(locTOFHitMatchParamsVector[loc_k]->dFlightTime);
501  tofList().setTflightvar(locTOFHitMatchParamsVector[loc_k]->dFlightTimeVariance);
502 
503  tofList().setDeltax(locTOFHitMatchParamsVector[loc_k]->dDeltaXToHit);
504  tofList().setDeltay(locTOFHitMatchParamsVector[loc_k]->dDeltaYToHit);
505  }
506 
507  vector<shared_ptr<const DSCHitMatchParams>> locSCHitMatchParamsVector;
508  locDetectorMatches[loc_i]->Get_SCMatchParams(tracks[loc_j], locSCHitMatchParamsVector);
509  for(size_t loc_k = 0; loc_k < locSCHitMatchParamsVector.size(); ++loc_k)
510  {
511  hddm_r::ScMatchParamsList scList = matches().addScMatchParamses(1);
512  scList().setTrack(loc_j);
513 
514  size_t locSCindex = 0;
515  for(; locSCindex < starthits.size(); ++locSCindex)
516  {
517  if(starthits[locSCindex] == locSCHitMatchParamsVector[loc_k]->dSCHit)
518  break;
519  }
520  scList().setHit(locSCindex);
521 
522  scList().setDEdx(locSCHitMatchParamsVector[loc_k]->dEdx);
523  scList().setDeltaphi(locSCHitMatchParamsVector[loc_k]->dDeltaPhiToHit);
524  scList().setEhit(locSCHitMatchParamsVector[loc_k]->dHitEnergy);
525  scList().setPathlength(locSCHitMatchParamsVector[loc_k]->dPathLength);
526  scList().setTflight(locSCHitMatchParamsVector[loc_k]->dFlightTime);
527  scList().setTflightvar(locSCHitMatchParamsVector[loc_k]->dFlightTimeVariance);
528  scList().setThit(locSCHitMatchParamsVector[loc_k]->dHitTime);
529  scList().setThitvar(locSCHitMatchParamsVector[loc_k]->dHitTimeVariance);
530  }
531 
532 
533  shared_ptr<const DDIRCMatchParams> locDIRCMatchParams;
534  map<shared_ptr<const DDIRCMatchParams>, vector<const DDIRCPmtHit*> > locDIRCTrackMatchParamsMap;
535  DDetectorMatches *locDetectorMatch = (DDetectorMatches*)locDetectorMatches[loc_i];
536  locDetectorMatch->Get_DIRCMatchParams(tracks[loc_j], locDIRCMatchParams);
537  locDetectorMatch->Get_DIRCTrackMatchParamsMap(locDIRCTrackMatchParamsMap);
538 
539  if(locDIRCMatchParams) {
540  hddm_r::DircMatchParamsList dircList = matches().addDircMatchParamses(1);
541  dircList().setTrack(loc_j);
542 
543  vector<const DDIRCPmtHit*> locDIRCHitTrackMatch = (vector<const DDIRCPmtHit*>)locDIRCTrackMatchParamsMap[locDIRCMatchParams];
544  for(size_t loc_k = 0; loc_k < locDIRCPmtHits.size(); ++loc_k) {
545  const DDIRCPmtHit* locDIRCPmtHit = (DDIRCPmtHit*)locDIRCPmtHits[loc_k];
546  if(find(locDIRCHitTrackMatch.begin(), locDIRCHitTrackMatch.end(), locDIRCPmtHit) != locDIRCHitTrackMatch.end()) {
547  hddm_r::DircMatchHitList dircHitList = matches().addDircMatchHits(1);
548  dircHitList().setTrack(loc_j);
549  dircHitList().setHit(loc_k);
550  }
551  }
552 
553  vector<DTrackFitter::Extrapolation_t> extrapolations=tracks[loc_j]->extrapolations.at(SYS_DIRC);
554  DVector3 locProjPos = locDIRCMatchParams->dExtrapolatedPos;
555  DVector3 locProjMom = locDIRCMatchParams->dExtrapolatedMom;
556  double locFlightTime = locDIRCMatchParams->dExtrapolatedTime;
557  dircList().setX(locProjPos.X());
558  dircList().setY(locProjPos.Y());
559  dircList().setZ(locProjPos.Z());
560  dircList().setT(locFlightTime);
561  dircList().setPx(locProjMom.X());
562  dircList().setPy(locProjMom.Y());
563  dircList().setPz(locProjMom.Z());
564  dircList().setExpectthetac(locDIRCMatchParams->dExpectedThetaC);
565  dircList().setThetac(locDIRCMatchParams->dThetaC);
566  dircList().setDeltat(locDIRCMatchParams->dDeltaT);
567  dircList().setLele(locDIRCMatchParams->dLikelihoodElectron);
568  dircList().setLpi(locDIRCMatchParams->dLikelihoodPion);
569  dircList().setLk(locDIRCMatchParams->dLikelihoodKaon);
570  dircList().setLp(locDIRCMatchParams->dLikelihoodProton);
571  dircList().setNphotons(locDIRCMatchParams->dNPhotons);
572  }
573 
574  double locFlightTimePCorrelation = 0.0;
575  if(locDetectorMatches[loc_i]->Get_FlightTimePCorrelation(tracks[loc_j], SYS_BCAL, locFlightTimePCorrelation))
576  {
577  hddm_r::TflightPCorrelationList correlationList = matches().addTflightPCorrelations(1);
578  correlationList().setTrack(loc_j);
579  correlationList().setSystem(SYS_BCAL);
580  correlationList().setCorrelation(locFlightTimePCorrelation);
581  }
582  if(locDetectorMatches[loc_i]->Get_FlightTimePCorrelation(tracks[loc_j], SYS_FCAL, locFlightTimePCorrelation))
583  {
584  hddm_r::TflightPCorrelationList correlationList = matches().addTflightPCorrelations(1);
585  correlationList().setTrack(loc_j);
586  correlationList().setSystem(SYS_FCAL);
587  correlationList().setCorrelation(locFlightTimePCorrelation);
588  }
589  if(locDetectorMatches[loc_i]->Get_FlightTimePCorrelation(tracks[loc_j], SYS_TOF, locFlightTimePCorrelation))
590  {
591  hddm_r::TflightPCorrelationList correlationList = matches().addTflightPCorrelations(1);
592  correlationList().setTrack(loc_j);
593  correlationList().setSystem(SYS_TOF);
594  correlationList().setCorrelation(locFlightTimePCorrelation);
595  }
596  if(locDetectorMatches[loc_i]->Get_FlightTimePCorrelation(tracks[loc_j], SYS_START, locFlightTimePCorrelation))
597  {
598  hddm_r::TflightPCorrelationList correlationList = matches().addTflightPCorrelations(1);
599  correlationList().setTrack(loc_j);
600  correlationList().setSystem(SYS_START);
601  correlationList().setCorrelation(locFlightTimePCorrelation);
602  }
603  }
604 
605  for(size_t loc_j = 0; loc_j < bcalshowers.size(); ++loc_j)
606  {
607  double locDeltaPhi = 0.0, locDeltaZ = 0.0;
608  if(!locDetectorMatches[loc_i]->Get_DistanceToNearestTrack(bcalshowers[loc_j], locDeltaPhi, locDeltaZ))
609  continue;
610 
611  hddm_r::BcalDOCAtoTrackList bcalDocaList = matches().addBcalDOCAtoTracks(1);
612  bcalDocaList().setShower(loc_j);
613  bcalDocaList().setDeltaphi(locDeltaPhi);
614  bcalDocaList().setDeltaz(locDeltaZ);
615  }
616 
617  for(size_t loc_j = 0; loc_j < fcalshowers.size(); ++loc_j)
618  {
619  double locDistance = 0.0;
620  if(!locDetectorMatches[loc_i]->Get_DistanceToNearestTrack(fcalshowers[loc_j], locDistance))
621  continue;
622 
623  hddm_r::FcalDOCAtoTrackList fcalDocaList = matches().addFcalDOCAtoTracks(1);
624  fcalDocaList().setShower(loc_j);
625  fcalDocaList().setDoca(locDistance);
626  }
627  }
628 
629  // write the resulting record to the output stream
630  bool locWriteStatus = Write_RESTEvent(locOutputFileName, locRecord);
631  locRecord.clear();
632  return locWriteStatus;
633 }
634 
635 string DEventWriterREST::Get_OutputFileName(string locOutputFileNameSubString) const
636 {
637  string locOutputFileName = dOutputFileBaseName;
638  if (locOutputFileNameSubString != "")
639  locOutputFileName += string("_") + locOutputFileNameSubString;
640  return (locOutputFileName + string(".hddm"));
641 }
642 
643 bool DEventWriterREST::Write_RESTEvent(string locOutputFileName, hddm_r::HDDM& locRecord) const
644 {
645  japp->WriteLock("RESTWriter");
646  {
647  //check to see if the REST file is open
648  if(Get_RESTOutputFilePointers().find(locOutputFileName) != Get_RESTOutputFilePointers().end())
649  {
650  //open: get pointer, write event
651  hddm_r::ostream* locOutputRESTFileStream = Get_RESTOutputFilePointers()[locOutputFileName].second;
652  japp->Unlock("RESTWriter");
653  *(locOutputRESTFileStream) << locRecord;
654  return true;
655  }
656 
657  //not open: open it
658  pair<ofstream*, hddm_r::ostream*> locRESTFilePointers(NULL, NULL);
659  locRESTFilePointers.first = new ofstream(locOutputFileName.c_str());
660  if(!locRESTFilePointers.first->is_open())
661  {
662  //failed to open
663  delete locRESTFilePointers.first;
664  japp->Unlock("RESTWriter");
665  return false;
666  }
667  locRESTFilePointers.second = new hddm_r::ostream(*locRESTFilePointers.first);
668 
669  // enable on-the-fly bzip2 compression on output stream
671  {
672  jout << " Enabling bz2 compression of output HDDM file stream" << std::endl;
673  locRESTFilePointers.second->setCompression(hddm_r::k_bz2_compression);
674  }
675  else
676  jout << " HDDM compression disabled" << std::endl;
677 
678  // enable a CRC data integrity check at the end of each event record
680  {
681  jout << " Enabling CRC data integrity check in output HDDM file stream" << std::endl;
682  locRESTFilePointers.second->setIntegrityChecks(hddm_r::k_crc32_integrity);
683  }
684  else
685  jout << " HDDM integrity checks disabled" << std::endl;
686 
687  // write a comment record at the head of the file
688  hddm_r::HDDM locCommentRecord;
689  hddm_r::ReconstructedPhysicsEventList res = locCommentRecord.addReconstructedPhysicsEvents(1);
690  hddm_r::CommentList comment = res().addComments();
691  comment().setText("This is a REST event stream...");
692  // write out any metadata if it's been set
693  if(HDDM_DATA_VERSION_STRING != "") {
694  hddm_r::DataVersionStringList dataVersionString = res().addDataVersionStrings();
695  dataVersionString().setText(HDDM_DATA_VERSION_STRING);
696  }
697  if(CCDB_CONTEXT_STRING != "") {
698  hddm_r::CcdbContextList ccdbContextString = res().addCcdbContexts();
699  ccdbContextString().setText(CCDB_CONTEXT_STRING);
700  }
701  *(locRESTFilePointers.second) << locCommentRecord;
702  locCommentRecord.clear();
703 
704  //write the event
705  *(locRESTFilePointers.second) << locRecord;
706 
707  //store the stream pointers
708  Get_RESTOutputFilePointers()[locOutputFileName] = locRESTFilePointers;
709  }
710  japp->Unlock("RESTWriter");
711 
712  return true;
713 }
714 
716 {
717  japp->WriteLock("RESTWriter");
718  {
720  if(Get_NumEventWriterThreads() > 0)
721  {
722  japp->Unlock("RESTWriter");
723  return; //not the last thread writing to REST files
724  }
725 
726  //last thread writing to REST files: close all files and free all memory
727  map<string, pair<ofstream*, hddm_r::ostream*> >::iterator locIterator;
728  for(locIterator = Get_RESTOutputFilePointers().begin(); locIterator != Get_RESTOutputFilePointers().end(); ++locIterator)
729  {
730  string locOutputFileName = locIterator->first;
731  if (locIterator->second.second != NULL)
732  delete locIterator->second.second;
733  if (locIterator->second.first != NULL)
734  delete locIterator->second.first;
735  std::cout << "Closed REST file " << locOutputFileName << std::endl;
736  }
737  Get_RESTOutputFilePointers().clear();
738  }
739  japp->Unlock("RESTWriter");
740 }
741 
742 int32_t DEventWriterREST::Convert_UnsignedIntToSigned(uint32_t locUnsignedInt) const
743 {
744  //Convert uint32_t to int32_t
745  //Scheme:
746  //If bit 32 is zero, then the int32_t is the same as the uint32_t: Positive or zero
747  //If bit 32 is one, and at least one other bit is 1, then the int32_t is -1 * uint32_t (after stripping the top bit)
748  //If bit 32 is one, and all other bits are zero, then the int32_t is the minimum int: -(2^31)
749  if((locUnsignedInt & 0x80000000) == 0)
750  return int32_t(locUnsignedInt); //bit 32 is zero: positive or zero
751 
752  //bit 32 is 1. see if there is another bit set
753  int32_t locTopBitStripped = int32_t(locUnsignedInt & uint32_t(0x7FFFFFFF)); //strip the top bit
754  if(locTopBitStripped == 0)
755  return numeric_limits<int32_t>::min(); //no other bit is set: minimum int
756  return -1*locTopBitStripped; //return the negative
757 }
DApplication * dapp
bool Get_DIRCTrackMatchParamsMap(map< shared_ptr< const DDIRCMatchParams >, vector< const DDIRCPmtHit * > > &locDIRCTrackMatchParamsMap)
TVector3 DVector3
Definition: DVector3.h:14
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
char string[256]
#define y
Definition: GlueX.h:29
Definition: GlueX.h:24
TLorentzVector DLorentzVector
static int PDGtype(Particle_t p)
Definition: GlueX.h:17
Definition: GlueX.h:19
JApplication * japp
bool Write_RESTEvent(JEventLoop *locEventLoop, string locOutputFileNameSubString) const
Definition: GlueX.h:20
DEventWriterREST(JEventLoop *locEventLoop, string locOutputFileBaseName)
Definition: GlueX.h:22
Definition: GlueX.h:26
int & Get_NumEventWriterThreads(void) const
map< string, pair< ofstream *, hddm_r::ostream * > > & Get_RESTOutputFilePointers(void) const
bool Get_DIRCMatchParams(const DTrackingData *locTrack, shared_ptr< const DDIRCMatchParams > &locMatchParams) const
TH1D * py[NCHANNELS]
static int matches(char *d, char *c, t_iostream *fp)
Definition: hddm_t.c:356
string HDDM_DATA_VERSION_STRING
int32_t Convert_UnsignedIntToSigned(uint32_t locUnsignedInt) const
string Get_OutputFileName(string locOutputFileNameSubString) const
Particle_t
Definition: particleType.h:12