StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StBTofNtupleMaker.cxx
1 /*******************************************************************
2  *
3  * $Id: StBTofNtupleMaker.cxx,v 1.3 2018/02/26 23:26:51 smirnovd Exp $
4  *
5  * Author: Xin Dong
6  *****************************************************************
7  *
8  * Description: example maker to get the matched TOFr cells and fill
9  * into TOFr tree.
10  *
11  *******************************************************************/
12 #include <iostream>
13 #include <math.h>
14 #include <vector>
15 #include <stdlib.h>
16 #include <stdio.h>
17 #include <algorithm>
18 #include <iterator>
19 
20 #include "StEventTypes.h"
21 #include "Stypes.h"
22 #include "StThreeVectorF.hh"
23 #include "StMeasuredPoint.h"
24 #include "StDedxPidTraits.h"
25 #include "StBTofPidTraits.h"
26 #include "StTrackPidTraits.h"
27 #include "StBTofPidTraits.h"
28 #include "StBTofCollection.h"
29 #include "StBTofHit.h"
30 #include "StBTofHeader.h"
31 #include "StarClassLibrary/StParticleTypes.hh"
32 #include "StarClassLibrary/StParticleDefinition.hh"
33 #include "StMuDSTMaker/COMMON/StMuUtilities.h"
34 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
35 #include "StMuDSTMaker/COMMON/StMuDst.h"
36 #include "StMuDSTMaker/COMMON/StMuEvent.h"
37 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
38 #include "StMuDSTMaker/COMMON/StMuTrack.h"
39 #include "StHelix.hh"
40 #include "StTrackGeometry.h"
41 #include "StDcaGeometry.h"
42 #include "StParticleTypes.hh"
43 #include "StTpcDedxPidAlgorithm.h"
44 #include "StHit.h"
45 #include "StAssociationMaker/StTrackPairInfo.hh"
46 #include "StEventUtilities/StuRefMult.hh"
47 #include "PhysicalConstants.h"
48 #include "StPhysicalHelixD.hh"
49 #include "TTree.h"
50 #include "TFile.h"
51 #include "TH1.h"
52 #include "TH2.h"
53 #include "StMemoryInfo.hh"
54 #include "StMessMgr.h"
55 #include "StTimer.hh"
56 #include "tables/St_g2t_vertex_Table.h" // tmp for Dz(vertex)
57 #include "tables/St_vertexSeed_Table.h" //
58 
59 #include "StBTofUtil/tofPathLength.hh"
60 #include "StBTofUtil/StBTofGeometry.h"
61 #include "StBTofUtil/StBTofDaqMap.h"
62 #include "StBTofUtil/StBTofHitCollection.h"
63 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
64 #include "StMuDSTMaker/COMMON/StMuDst.h"
65 #include "StMuDSTMaker/COMMON/StMuBTofHit.h"
66 #include "StMuDSTMaker/COMMON/StMuTrack.h"
67 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
68 #include "StMuDSTMaker/COMMON/StMuBTofPidTraits.h"
69 #include "StMuDSTMaker/COMMON/StMuTriggerIdCollection.h"
70 
71 
72 #include "StBTofNtupleMaker.h"
73 
74 #include "StEnumerations.h"
75 
76 
77 //---------------------------------------------------------------------------
79 StBTofNtupleMaker::StBTofNtupleMaker(const Char_t *name="tofNtuple", const Char_t *outname="tofntuple.root") : StMaker(name) {
80  mTupleFileName=outname;
81 
82  doPrintMemoryInfo = kFALSE;
83  doPrintCpuInfo = kFALSE;
84 
85  //setOuterGeometry(kFALSE);
86  mBeamHelix = 0;
87 }
88 
90 StBTofNtupleMaker::~StBTofNtupleMaker(){ /* nope */}
91 
92 //---------------------------------------------------------------------------
94 Int_t StBTofNtupleMaker::Init(){
95 
96  if(!mUseEventVertex) { // if not using, use beamHelix to caculate dca
97  StThreeVectorD MomFstPt(0.,0.,9999.);
98  StThreeVectorD origin(0.,0.,0.);
99  mBeamHelix = new StPhysicalHelixD(MomFstPt,origin,0.5*tesla,1.);
100  }
101 
102  if (mTupleFileName!="") bookNtuples();
103 
104  mAcceptedEvents = 0;
105  mPvpdEntries = 0;
106  mBTofEvents = 0;
107  mBTofEntries = 0;
108 
109  return kStOK;
110 }
111 
112 
113 Int_t StBTofNtupleMaker::InitRun(int runnumber) {
114 
115  if(mUseEventVertex) return kStOK;
116 
117  //========== Set Beam Line =====================
118  double x0 = 0.;
119  double y0 = 0.;
120  double dxdz = 0.;
121  double dydz = 0.;
122 
123  //Get Current Beam Line Constraint from database
124  TDataSet* dbDataSet = this->GetDataBase("Calibrations/rhic");
125 
126  if (dbDataSet) {
127  vertexSeed_st* vSeed = ((St_vertexSeed*) (dbDataSet->FindObject("vertexSeed")))->GetTable();
128 
129  x0 = vSeed->x0;
130  y0 = vSeed->y0;
131  dxdz = vSeed->dxdz;
132  dydz = vSeed->dydz;
133  }
134  else {
135  LOG_INFO << "StTofrNtupleMaker -- No Database for beamline" << endm;
136  }
137 
138  LOG_INFO << "BeamLine Constraint (StTofrNtupleMaker): " << endm;
139  LOG_INFO << "x(z) = " << x0 << " + " << dxdz << " * z" << endm;
140  LOG_INFO << "y(z) = " << y0 << " + " << dydz << " * z" << endm;
141 
142  //set by hand
143 // x0=0.;
144 // y0=0.;
145  StThreeVectorD origin(x0,y0,0.0);
146  double pt = 88889999;
147  double nxy=::sqrt(dxdz*dxdz + dydz*dydz);
148  if(nxy<1.e-5){ // beam line _MUST_ be tilted
149  LOG_WARN << "StTofrNtupleMaker:: Beam line must be tilted!" << endm;
150  nxy=dxdz=1.e-5;
151  }
152  double p0=pt/nxy;
153  double px = p0*dxdz;
154  double py = p0*dydz;
155  double pz = p0; // approximation: nx,ny<<0
156  StThreeVectorD MomFstPt(px*GeV, py*GeV, pz*GeV);
157  //delete mBeamHelix;
158  mBeamHelix = new StPhysicalHelixD(MomFstPt,origin,0.5*tesla,1.);
159 
160  return kStOK;
161 }
162 
163 Int_t StBTofNtupleMaker::FinishRun(int runnumber)
164 {
165  if(mBeamHelix) delete mBeamHelix;
166  mBeamHelix = 0;
167 
168  return kStOK;
169 }
170 
172 Int_t StBTofNtupleMaker::Finish() {
173 
174  if (!(mTupleFileName=="")){
175  mTupleFile->Write();
176  mTupleFile->Close();
177  LOG_INFO << "StBTofNtupleMaker::Finish() ntuple file "
178  << mTupleFileName << " closed." << endm;
179  }
180 
181  //delete mPvpdTuple;
182  //delete mCellTuple;
183  //delete mTupleFile;
184 
185  LOG_INFO << "StBTofNtupleMaker -- statistics" << endm;
186  LOG_INFO << " accepted events : " << mAcceptedEvents << endm;
187  LOG_INFO << " pVPD entries : " << mPvpdEntries << endm;
188  LOG_INFO << " BTof entries/events : " << mBTofEntries << "/" << mBTofEvents << endm;
189  return kStOK;
190 }
191 
192 
193 //---------------------------------------------------------------------------
195 Int_t StBTofNtupleMaker::Make(){
196  LOG_INFO << "StBTofNtupleMaker -- welcome" << endm;
197 
198  if(!mMuDstIn) processStEvent();
199  else processMuDst();
200 
201  return kStOK;
202 }
203 
204 //---------------------------------------------------------------------------
205 void StBTofNtupleMaker::processStEvent() {
206 
207  StEvent *mEvent = (StEvent *) GetInputDS("StEvent");
208  //.........................................................................
209  // event selection ...
210  if (!mEvent||!mEvent->btofCollection()) {
211  LOG_INFO << "StBTofNtupleMaker -- nothing to do ... bye-bye"<< endm;
212  return;
213  }
214 
215  mAcceptedEvents++;
216  StTimer timer;
217  if (doPrintCpuInfo) timer.start();
218  if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
219 
220  //.........................................................................
221  // Collect global data for both ntuples
222 
223  StEventInfo *info = mEvent->info();
224  if(info) {
225  if(Debug()) { LOG_INFO<<"runId: "<<mEvent->runId()<<" evtId: "<<mEvent->id()<<endm; }
226  }
227 
228  StThreeVectorD pVtx(-999., -999., -999.);
229  if(mEvent->primaryVertex()) {
230  pVtx = mEvent->primaryVertex()->position();
231  }
232  //mCellData.trgword = triggerWord;
233 
234  mCellData.vertexX = pVtx.x();
235  mCellData.vertexY = pVtx.y();
236  mCellData.vertexZ = pVtx.z();
237 
238  //-------- fill mEvent summary info -----------
239  if(info) {
240  mCellData.run = mEvent->runId(); // the run number
241  mCellData.evt = mEvent->id(); // the event number
242  } else {
243  mCellData.run = 0;
244  mCellData.evt = 0;
245  }
246 
247  //-- read in TOF info
248  int ntofhits = 0;
249  StBTofCollection *theTof = mEvent->btofCollection();
250  if(Debug()&&theTof) { LOG_INFO << "got btof Collection"<<endm; }
251 
252  StBTofHeader* tofHeader = theTof->tofHeader();
253  if(!tofHeader) {
254  LOG_WARN << " No TOF header ... bye-bye" << endm;
255  return;
256  }
257  if(Debug()&&tofHeader) { LOG_INFO << "got tof Header"<<endm; }
258  mCellData.tStart = tofHeader->tStart();
259  mCellData.tDiff = tofHeader->tDiff();
260  mCellData.vpdVz = tofHeader->vpdVz();
261 
262  //initialize vpd content
263  for(int i=0;i<19;i++){
264  mCellData.vpdLeEast[i] = 0;
265  mCellData.vpdTotEast[i] = 0;
266  mCellData.vpdLeWest[i] = 0;
267  mCellData.vpdTotWest[i] = 0;
268  }
269 
270  unsigned int vpdEast=0, vpdWest=0, nVpdEast=0, nVpdWest=0;
271 
272  if (theTof && theTof->hitsPresent()){
273 
274  StSPtrVecBTofHit& hits = theTof->tofHits();
275  if(Debug()) { LOG_INFO << hits.size() << " hits"<<endm; }
276  for(size_t i=0;i<hits.size();i++) {
277  StBTofHit *aHit = (StBTofHit *)hits[i];
278  int trayId = aHit->tray();
279  StThreeVector<double> globalPos;
280  if(Debug()) { LOG_INFO << "tray Id = "<<trayId<<endm; }
281  if(trayId==122){//vpd East
282  int tubeId = aHit->cell()-1;
283  mCellData.vpdLeEast[tubeId] = aHit->leadingEdgeTime();
284  mCellData.vpdTotEast[tubeId] = aHit->tot();
285  vpdEast += 1<<tubeId;
286  nVpdEast++;
287 
288  } else if(trayId==121) {//vpd West
289  int tubeId = aHit->cell()-1;
290  mCellData.vpdLeWest[tubeId] = aHit->leadingEdgeTime();
291  mCellData.vpdTotWest[tubeId] = aHit->tot();
292  vpdWest += 1<<tubeId;
293  nVpdWest++;
294  }
295  else if(trayId<=120&&trayId>=0) {//TOF
296  mCellData.tray[ntofhits] = aHit->tray();
297  mCellData.module[ntofhits] = aHit->module();
298  mCellData.cell[ntofhits] = aHit->cell();
299  mCellData.leTime[ntofhits] = aHit->leadingEdgeTime();
300  mCellData.tot[ntofhits] = aHit->tot();
301 
302  if(Debug()) { LOG_INFO <<"tray/module/cell/letime/tot="<<aHit->tray()<<"/"<<aHit->module()<<"/"<<aHit->cell()<<"/"<<aHit->leadingEdgeTime()<<"/"<<aHit->tot()<<endm; }
303  //- track information
304  StTrack *thisTrack = aHit->associatedTrack();
305  if(!thisTrack) continue;
306  if(Debug()&&thisTrack) { LOG_INFO <<" got associated track"<<endm; }
307  StGlobalTrack* globalTrack = (StGlobalTrack*)thisTrack->node()->track(global);
308  if(Debug()) { LOG_INFO << "got global track"<<endm; }
309 
310  StTrackGeometry *theDcaGeometry = 0;
311  if(mUseEventVertex) {
312  StPrimaryTrack* pTrack = (StPrimaryTrack*)thisTrack->node()->track(primary);
313  if(!pTrack) continue; // save only primary tracks
314  if(pTrack->vertex()!=mEvent->primaryVertex()) continue; // first primary vertex
315  theDcaGeometry = pTrack->geometry();
316  } else {
317  theDcaGeometry = globalTrack->geometry();
318  }
319 
320  if(Debug()) { LOG_INFO << "got DcaGeometry"<<endm; }
321  StThreeVectorF momentum = theDcaGeometry->momentum();
322 
323  mCellData.trackId[ntofhits] = (Int_t) thisTrack->key();
324  mCellData.charge[ntofhits] = theDcaGeometry->charge();
325  mCellData.pt[ntofhits] = momentum.perp();
326  mCellData.eta[ntofhits] = momentum.pseudoRapidity();
327  mCellData.phi[ntofhits] = momentum.phi();
328 
329 
330  //-- dEdx and Tof PID traits
331  StSPtrVecTrackPidTraits& traits = thisTrack->pidTraits();
332  for (unsigned int it=0;it<traits.size();it++){
333  if (traits[it]->detector() == kTpcId){
334  StDedxPidTraits *pid = dynamic_cast<StDedxPidTraits*>(traits[it]);
335  if (pid && pid->method() ==kTruncatedMeanId){
336  mCellData.dedx[ntofhits] = pid->mean()*1e6;
337  mCellData.nHitsDedx[ntofhits] = pid->numberOfPoints();
338  }
339  } else if (traits[it]->detector() == kTofId) {
340  StBTofPidTraits* tofpid = dynamic_cast<StBTofPidTraits*>(traits[it]);
341  if(tofpid){
342  mCellData.matchFlag[ntofhits] = tofpid->matchFlag();
343  mCellData.yLocal[ntofhits] = tofpid->yLocal();
344  mCellData.zLocal[ntofhits] = tofpid->zLocal();
345  mCellData.thetaLocal[ntofhits] = tofpid->thetaLocal();
346  globalPos = tofpid->position();
347  mCellData.xGlobal[ntofhits] = globalPos.x();
348  mCellData.yGlobal[ntofhits] = globalPos.y();
349  mCellData.zGlobal[ntofhits] = globalPos.z();
350 
351  mCellData.tofCorr[ntofhits] = tofpid->timeOfFlight();
352  mCellData.length[ntofhits] = tofpid->pathLength();
353  mCellData.beta[ntofhits] = tofpid->beta();
354  }
355  }
356  }
357 
358  if(thisTrack->detectorInfo()) {
359  mCellData.nHits[ntofhits] = thisTrack->detectorInfo()->numberOfPoints(kTpcId);
360  } else {
361  mCellData.nHits[ntofhits] = 0;
362  }
363  mCellData.nHitsFit[ntofhits] = thisTrack->fitTraits().numberOfFitPoints(kTpcId);
364 
365 
366  //-- nSigma
367  static StTpcDedxPidAlgorithm PidAlgorithm;
368  static StElectron* Electron = StElectron::instance();
369  static StPionPlus* Pion = StPionPlus::instance();
370  static StKaonPlus* Kaon = StKaonPlus::instance();
371  static StProton* Proton = StProton::instance();
372  const StParticleDefinition* pd = thisTrack->pidTraits(PidAlgorithm);
373 
374  if (pd) {
375  mCellData.nSigE[ntofhits] = fabsMin(PidAlgorithm.numberOfSigma(Electron), __SIGMA_SCALE__);
376  mCellData.nSigPi[ntofhits] = fabsMin(PidAlgorithm.numberOfSigma(Pion),__SIGMA_SCALE__);
377  mCellData.nSigK[ntofhits] = fabsMin(PidAlgorithm.numberOfSigma(Kaon),__SIGMA_SCALE__);
378  mCellData.nSigP[ntofhits] = fabsMin(PidAlgorithm.numberOfSigma(Proton),__SIGMA_SCALE__);
379  }
380 
381  //-- project track onto beam line
382  if(mUseEventVertex) {
383  double s = globalTrack->geometry()->helix().pathLength(pVtx, false);
384  StThreeVectorD dca3D = globalTrack->geometry()->helix().at(s) - pVtx;
385  mCellData.dcaX[ntofhits] = dca3D.x();
386  mCellData.dcaY[ntofhits] = dca3D.y();
387  mCellData.dcaZ[ntofhits] = dca3D.z();
388  //-- get path length
389  mCellData.length[ntofhits] = tofPathLength(&pVtx, &globalPos, theDcaGeometry->helix().curvature());
390  } else {
391  StThreeVector<double> tofPos = theDcaGeometry->helix().at(theDcaGeometry->helix().pathLengths(*mBeamHelix).first);
392  StThreeVector<double> beamPos = mBeamHelix->at(theDcaGeometry->helix().pathLengths(*mBeamHelix).second);
393  StThreeVector<double> dcatof = tofPos - beamPos;
394  mCellData.dcaX[ntofhits] = dcatof.x();
395  mCellData.dcaY[ntofhits] = dcatof.y();
396  mCellData.dcaZ[ntofhits] = tofPos.z();
397  if(Debug()) {
398  LOG_INFO<<" tofPos(x,y,z) = "<<tofPos.x()<<","<<tofPos.y()<<","<<tofPos.z()<<endm;
399  LOG_INFO<<" beamPos(x,y,z) = "<<beamPos.x()<<","<<beamPos.y()<<","<<beamPos.z()<<endm;
400  LOG_INFO<<" dca (x,y,z) = "<<dcatof.x()<<","<<dcatof.y()<<","<<dcatof.z()<<endm;
401  LOG_INFO<<" 2D dca = "<<sqrt(pow(dcatof.x(),2)+pow(dcatof.y(),2))<<endm;
402  LOG_INFO<<" 2D signed dca = "<<theDcaGeometry->helix().geometricSignedDistance(beamPos.x(),beamPos.y())<<endm;
403  }
404  //-- get path length
405  mCellData.length[ntofhits] = tofPathLength(&tofPos, &globalPos, theDcaGeometry->helix().curvature());
406  }
407 
408  if(Debug()) { LOG_INFO << "pathLength="<< mCellData.length[ntofhits] <<endm; }
409 
410  ntofhits++;
411  }
412  }
413  }
414  mCellData.nTofHits = ntofhits;
415  mCellData.vpdEast = vpdEast;
416  mCellData.vpdWest = vpdWest;
417  mCellData.numberOfVpdEast = nVpdEast;
418  mCellData.numberOfVpdWest = nVpdWest;
419 
420  if(Debug()) { LOG_INFO << " Three are " << ntofhits << " tof hits in this event! " << endm; }
421  mBTofEntries = ntofhits;
422  mCellTuple->Fill();
423 
424  //- debug info`
425  if (doPrintMemoryInfo) {
426  StMemoryInfo::instance()->snapshot();
427  StMemoryInfo::instance()->print();
428  }
429  if (doPrintCpuInfo) {
430  timer.stop();
431  LOG_INFO << "CPU time for StEventMaker::Make(): "
432  << timer.elapsedTime() << " sec\n" << endm;
433  }
434 }
435 
436 //---------------------------------------------------------------------------
437 void StBTofNtupleMaker::processMuDst() {
438 
439  StMuDstMaker *mudstMaker = (StMuDstMaker*) GetMaker("MuDst");
440  if(!mudstMaker) {
441  LOG_WARN << " No MuDstMaker ... bye-bye" << endm;
442  return;
443  }
444  mMuDst = mudstMaker->muDst();
445  if(!mMuDst) {
446  LOG_WARN << " No MuDst ... bye-bye" << endm;
447  return;
448  }
449 
450  mAcceptedEvents++;
451  StTimer timer;
452  if (doPrintCpuInfo) timer.start();
453  if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
454 
455  //.........................................................................
456  // Collect global data for both ntuples
457 
458  StMuEvent *mMuEvent = mMuDst->event();
459  Bool_t istrigger = 0;
460 
461 
462  if(mMuEvent) {
463  if(Debug()) { LOG_INFO<<"runId: "<<mMuEvent->runId()<<" evtId: "<<mMuEvent->eventId()<<endm; }
464  mCellData.run = mMuEvent->runId(); // the run number
465  mCellData.evt = mMuEvent->eventId(); // the event number
466  istrigger = ((mMuEvent->triggerIdCollection().nominal().isTrigger(7))||(mMuEvent->triggerIdCollection().nominal().isTrigger(250107)));
467  } else {
468  mCellData.run = 0;
469  mCellData.evt = 0;
470  }
471 // if(istrigger<1) return;
472 
473  StThreeVectorD pVtx(-999., -999., -999.);
474  if(mMuDst->primaryVertex()) {
475  pVtx = mMuDst->primaryVertex()->position();
476  }
477  mCellData.vertexX = pVtx.x();
478  mCellData.vertexY = pVtx.y();
479  mCellData.vertexZ = pVtx.z();
480 
481  //-- read in TOF info
482  StBTofHeader* tofHeader = mMuDst->btofHeader();
483  if(!tofHeader) {
484  LOG_WARN << " No TOF Header ... bye-bye" << endm;
485  return;
486  }
487  if(Debug()&&tofHeader) { LOG_INFO << "got tof Header"<<endm; }
488  mCellData.tStart = tofHeader->tStart();
489  mCellData.tDiff = tofHeader->tDiff();
490  mCellData.vpdVz = tofHeader->vpdVz();
491 
492  //initialize vpd content
493  for(int i=0;i<19;i++){
494  mCellData.vpdLeEast[i] = 0;
495  mCellData.vpdTotEast[i] = 0;
496  mCellData.vpdLeWest[i] = 0;
497  mCellData.vpdTotWest[i] = 0;
498  }
499 
500  unsigned int vpdEast=0, vpdWest=0, nVpdEast=0, nVpdWest=0;
501  int nMax = mMuDst->numberOfBTofHit();
502  int ntofhits = 0;
503  if(Debug()) { LOG_INFO << nMax << " hits"<<endm; }
504  for(int i=0;i<nMax;i++) {
505  StMuBTofHit *aHit = (StMuBTofHit *)mMuDst->btofHit(i);
506  int trayId = aHit->tray();
507  StThreeVector<double> globalPos;
508  if(Debug()) { LOG_INFO << "tray Id = "<<trayId<<endm; }
509  if(trayId==122){//vpd East
510  int tubeId = aHit->cell()-1;
511  mCellData.vpdLeEast[tubeId] = aHit->leadingEdgeTime();
512  mCellData.vpdTotEast[tubeId] = aHit->tot();
513  vpdEast += 1<<tubeId;
514  nVpdEast++;
515 
516  } else if(trayId==121) {//vpd West
517  int tubeId = aHit->cell()-1;
518  mCellData.vpdLeWest[tubeId] = aHit->leadingEdgeTime();
519  mCellData.vpdTotWest[tubeId] = aHit->tot();
520  vpdWest += 1<<tubeId;
521  nVpdWest++;
522  }
523  else if(trayId<=120&&trayId>=0) {//TOF
524  StMuTrack *globalTrack = aHit->globalTrack();
525  if(!globalTrack) continue;
526  if(Debug()) { LOG_INFO << "got global track"<<endm; }
527 
528  mCellData.tray[ntofhits] = aHit->tray();
529  mCellData.module[ntofhits] = aHit->module();
530  mCellData.cell[ntofhits] = aHit->cell();
531  mCellData.leTime[ntofhits] = aHit->leadingEdgeTime();
532  mCellData.tot[ntofhits] = aHit->tot();
533 
534  if(Debug()) { LOG_INFO <<"tray/module/cell/letime/tot="<<aHit->tray()<<"/"<<aHit->module()<<"/"<<aHit->cell()<<"/"<<aHit->leadingEdgeTime()<<"/"<<aHit->tot()<<endm; }
535 
536  StThreeVectorF momentum;
537  StMuTrack *pTrack = 0;
538  if(mUseEventVertex) {
539  pTrack = aHit->primaryTrack();
540  if(!pTrack) continue;
541  if(pTrack->vertexIndex()!=0) continue; // only select the first one
542  momentum = pTrack->momentum();
543  } else {
544  momentum = globalTrack->momentum();
545  }
546 
547  mCellData.trackId[ntofhits] = (Int_t) globalTrack->id();
548  mCellData.charge[ntofhits] = globalTrack->charge();
549  mCellData.pt[ntofhits] = momentum.perp();
550  mCellData.eta[ntofhits] = momentum.pseudoRapidity();
551  mCellData.phi[ntofhits] = momentum.phi();
552 
553  mCellData.dedx[ntofhits] = globalTrack->dEdx()*1e6;
554  mCellData.nHitsDedx[ntofhits] = globalTrack->nHitsDedx();
555 
556  StMuBTofPidTraits tofpid = globalTrack->btofPidTraits();
557 
558  mCellData.matchFlag[ntofhits] = tofpid.matchFlag();
559  mCellData.yLocal[ntofhits] = tofpid.yLocal();
560  mCellData.zLocal[ntofhits] = tofpid.zLocal();
561  mCellData.thetaLocal[ntofhits] = tofpid.thetaLocal();
562  globalPos = tofpid.position();
563  mCellData.xGlobal[ntofhits] = globalPos.x();
564  mCellData.yGlobal[ntofhits] = globalPos.y();
565  mCellData.zGlobal[ntofhits] = globalPos.z();
566 
567  mCellData.tofCorr[ntofhits] = tofpid.timeOfFlight();
568  mCellData.length[ntofhits] = tofpid.pathLength();
569  mCellData.beta[ntofhits] = tofpid.beta();
570 
571  mCellData.nHits[ntofhits] = globalTrack->nHits();
572  mCellData.nHitsFit[ntofhits] = globalTrack->nHitsFit(kTpcId);
573 
574  mCellData.nSigE[ntofhits] = globalTrack->nSigmaElectron();
575  mCellData.nSigPi[ntofhits] = globalTrack->nSigmaPion();
576  mCellData.nSigK[ntofhits] = globalTrack->nSigmaKaon();
577  mCellData.nSigP[ntofhits] = globalTrack->nSigmaProton();
578 
579  //-- project track onto beam line
580  if(mUseEventVertex) {
581  mCellData.dcaX[ntofhits] = pTrack->dcaGlobal().x();
582  mCellData.dcaY[ntofhits] = pTrack->dcaGlobal().y();
583  mCellData.dcaZ[ntofhits] = pTrack->dcaGlobal().z();
584  //-- get path length
585  mCellData.length[ntofhits] = tofPathLength(&pVtx, &globalPos, pTrack->helix().curvature());
586  } else {
587  StPhysicalHelixD helix = globalTrack->helix();
588  StThreeVector<double> tofPos = helix.at(helix.pathLengths(*mBeamHelix).first);
589  StThreeVector<double> beamPos = mBeamHelix->at(helix.pathLengths(*mBeamHelix).second);
590  StThreeVector<double> dcatof = tofPos - beamPos;
591  mCellData.dcaX[ntofhits] = dcatof.x();
592  mCellData.dcaY[ntofhits] = dcatof.y();
593  mCellData.dcaZ[ntofhits] = tofPos.z();
594  if(Debug()) {
595  LOG_INFO<<" tofPos(x,y,z) = "<<tofPos.x()<<","<<tofPos.y()<<","<<tofPos.z()<<endm;
596  LOG_INFO<<"beamPos(x,y,z) = "<<beamPos.x()<<","<<beamPos.y()<<","<<beamPos.z()<<endm;
597  LOG_INFO<<" dca (x,y,z) = "<<dcatof.x()<<","<<dcatof.y()<<","<<dcatof.z()<<endm;
598  LOG_INFO<<" 2D dca = "<<sqrt(pow(dcatof.x(),2)+pow(dcatof.y(),2))<<endm;
599  LOG_INFO<<" 2D signed dca = "<<helix.geometricSignedDistance(beamPos.x(),beamPos.y())<<endm;
600  }
601  //-- get path length
602  mCellData.length[ntofhits] = tofPathLength(&tofPos, &globalPos, helix.curvature());
603  }
604 
605  if(Debug()) LOG_INFO << "pathLength="<< mCellData.length[ntofhits] <<endm;
606 
607  ntofhits++;
608  }
609  }
610  mCellData.nTofHits = ntofhits;
611  mCellData.vpdEast = vpdEast;
612  mCellData.vpdWest = vpdWest;
613  mCellData.numberOfVpdEast = nVpdEast;
614  mCellData.numberOfVpdWest = nVpdWest;
615 
616  if(Debug()) { LOG_INFO << " Three are " << ntofhits << " tof hits in this event! " << endm; }
617  mBTofEntries = ntofhits;
618  mCellTuple->Fill();
619 
620  //- debug info`
621  if (doPrintMemoryInfo) {
622  StMemoryInfo::instance()->snapshot();
623  StMemoryInfo::instance()->print();
624  }
625  if (doPrintCpuInfo) {
626  timer.stop();
627  LOG_INFO << "CPU time for StEventMaker::Make(): "
628  << timer.elapsedTime() << " sec\n" << endm;
629  }
630 }
631 
632 //---------------------------------------------------------------------------
634 void StBTofNtupleMaker::bookNtuples(){
635  mTupleFile = new TFile(mTupleFileName.c_str(), "RECREATE");
636  LOG_INFO << "StBTofNtupleMaker::bookNtuples() file "
637  << mTupleFileName << " opened" << endm;
638 
639  // BTof calibration ntuple
640  mCellTuple = new TTree("tof","BTof cell data");
641  mCellTuple->SetAutoSave(1000);
642  mCellTuple->Branch("run",&mCellData.run,"run/I");
643  mCellTuple->Branch("evt",&mCellData.evt,"evt/I");
644  // mCellTuple->Branch("trgword",&mCellData.trgword,"trgword/I");
645  mCellTuple->Branch("vertexX",&mCellData.vertexX,"vertexX/F");
646  mCellTuple->Branch("vertexY",&mCellData.vertexY,"vertexY/F");
647  mCellTuple->Branch("vertexZ",&mCellData.vertexZ,"vertexZ/F");
648  mCellTuple->Branch("vpdEast",&mCellData.vpdEast,"vpdEast/I");
649  mCellTuple->Branch("vpdWest",&mCellData.vpdWest,"vpdWest/I");
650  mCellTuple->Branch("numberOfVpdEast",&mCellData.numberOfVpdEast,"numberOfVpdEast/I");
651  mCellTuple->Branch("numberOfVpdWest",&mCellData.numberOfVpdWest,"numberOfVpdWest/I");
652  mCellTuple->Branch("tDiff",&mCellData.tDiff,"tDiff/F");
653  mCellTuple->Branch("tStart",&mCellData.tStart,"tStart/D");
654  mCellTuple->Branch("vpdVz",&mCellData.vpdVz,"vpdVz/F");
655 
656  mCellTuple->Branch("vpdLeEast",&mCellData.vpdLeEast,"vpdLeEast[19]/D");
657  mCellTuple->Branch("vpdLeWest",&mCellData.vpdLeWest,"vpdLeWest[19]/D");
658  mCellTuple->Branch("vpdTotEast",&mCellData.vpdTotEast,"vpdTotEast[19]/D");
659  mCellTuple->Branch("vpdTotWest",&mCellData.vpdTotWest,"vpdTotWest[19]/D");
660  mCellTuple->Branch("nTofHits",&mCellData.nTofHits,"nTofHits/I");
661  mCellTuple->Branch("tray",&mCellData.tray,"tray[nTofHits]/I");
662  mCellTuple->Branch("module",&mCellData.module,"module[nTofHits]/I");
663  mCellTuple->Branch("cell",&mCellData.cell,"cell[nTofHits]/I");
664  mCellTuple->Branch("leTime",&mCellData.leTime,"leTime[nTofHits]/D");
665  mCellTuple->Branch("tot",&mCellData.tot,"tot[nTofHits]/D");
666  mCellTuple->Branch("matchFlag",&mCellData.matchFlag,"matchFlag/I");
667  mCellTuple->Branch("yLocal",&mCellData.yLocal,"yLocal[nTofHits]/F");
668  mCellTuple->Branch("zLocal",&mCellData.zLocal,"zLocal[nTofHits]/F");
669  mCellTuple->Branch("thetaLocal",&mCellData.thetaLocal,"thetaLocal[nTofHits]/F");
670  mCellTuple->Branch("xGlobal",&mCellData.xGlobal,"xGlobal[nTofHits]/F");
671  mCellTuple->Branch("yGlobal",&mCellData.yGlobal,"yGlobal[nTofHits]/F");
672  mCellTuple->Branch("zGlobal",&mCellData.zGlobal,"zGlobal[nTofHits]/F");
673  mCellTuple->Branch("trackId",&mCellData.trackId,"trackId[nTofHits]/I");
674  mCellTuple->Branch("charge",&mCellData.charge,"charge[nTofHits]/I");
675  mCellTuple->Branch("pt",&mCellData.pt,"pt[nTofHits]/F");
676  mCellTuple->Branch("eta",&mCellData.eta,"eta[nTofHits]/F");
677  mCellTuple->Branch("phi",&mCellData.phi,"phi[nTofHits]/F");
678  mCellTuple->Branch("dcaX",&mCellData.dcaX,"dcaX[nTofHits]/F");
679  mCellTuple->Branch("dcaY",&mCellData.dcaY,"dcaY[nTofHits]/F");
680  mCellTuple->Branch("dcaZ",&mCellData.dcaZ,"dcaZ[nTofHits]/F");
681  mCellTuple->Branch("length",&mCellData.length,"length[nTofHits]/F");
682  mCellTuple->Branch("nHits",&mCellData.nHits,"nHits[nTofHits]/I");
683  mCellTuple->Branch("nHitsFit",&mCellData.nHitsFit,"nHitsFit[nTofHits]/I");
684  mCellTuple->Branch("nHitsDedx",&mCellData.nHitsDedx,"nHitsDedx[nTofHits]/I");
685  mCellTuple->Branch("dedx",&mCellData.dedx,"dedx[nTofHits]/F");
686  mCellTuple->Branch("nSigE",&mCellData.nSigE,"nSigE[nTofHits]/F");
687  mCellTuple->Branch("nSigPi",&mCellData.nSigPi,"nSigPi[nTofHits]/F");
688  mCellTuple->Branch("nSigK",&mCellData.nSigK,"nSigK[nTofHits]/F");
689  mCellTuple->Branch("nSigP",&mCellData.nSigP,"nSigP[nTofHits]/F");
690  mCellTuple->Branch("tofCorr",&mCellData.tofCorr,"tofCorr[nTofHits]/F");
691  mCellTuple->Branch("beta",&mCellData.beta,"beta[nTofHits]/F");
692 
693  return;
694 }
695 
696 
697 /*****************************************************************
698  *
699  * $Log: StBTofNtupleMaker.cxx,v $
700  * Revision 1.3 2018/02/26 23:26:51 smirnovd
701  * StTof: Remove outdated ClassImp macro
702  *
703  * Revision 1.2 2018/02/26 23:13:20 smirnovd
704  * Move embedded CVS log messages to the end of file
705  *
706  * Revision 1.1 2010/04/09 00:28:48 dongx
707  * First release
708  *
709  * Revision 1.1 2010/04/09 00:16:05 dongx
710  * first release
711  *
712  */
StMuDst::primaryVertex
static StMuPrimaryVertex * primaryVertex()
return pointer to current primary vertex
Definition: StMuDst.h:322
StMuDst::btofHit
static StMuBTofHit * btofHit(int i)
returns pointer to the i-th muBTofHit
Definition: StMuDst.h:419
StMuTrack::vertexIndex
Int_t vertexIndex() const
Returns index of associated primary vertex.
Definition: StMuTrack.cxx:600
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
StHelix::pathLengths
pair< double, double > pathLengths(const StHelix &, double minStepSize=10 *micrometer, double minRange=10 *centimeter) const
path lengths at dca between two helices
Definition: StHelix.cc:511
StMuTrack::id
short id() const
Returns the track id(or key), is unique for a track node, i.e. global and primary tracks have the sam...
Definition: StMuTrack.h:228
StMuTrack::nHitsDedx
UShort_t nHitsDedx() const
Return number of hits used for dEdx.
Definition: StMuTrack.h:238
StMuBTofPidTraits::matchFlag
unsigned char matchFlag() const
Matching information.
Definition: StMuBTofPidTraits.h:115
StBTofNtupleMaker::~StBTofNtupleMaker
~StBTofNtupleMaker()
default empty destructor
Definition: StBTofNtupleMaker.cxx:90
StHelix::pathLength
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
StMuTrack::helix
StPhysicalHelixD helix() const
Returns inner helix (first measured point)
Definition: StMuTrack.cxx:407
StMuTrack::nHitsFit
UShort_t nHitsFit() const
Return total number of hits used in fit.
Definition: StMuTrack.h:239
StMuTrack::charge
Short_t charge() const
Returns charge.
Definition: StMuTrack.h:255
StBTofPidTraits::matchFlag
unsigned char matchFlag() const
Matching information.
Definition: StBTofPidTraits.h:120
StMuTrack::nSigmaPion
Double_t nSigmaPion() const
Returns Craig&#39;s distance to the calculated dE/dx band for pions in units of sigma.
Definition: StMuTrack.h:245
StMuTrack::nSigmaElectron
Double_t nSigmaElectron() const
Returns Craig&#39;s distance to the calculated dE/dx band for electrons in units of sigma.
Definition: StMuTrack.h:244
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
StMuTrack::momentum
const StThreeVectorF & momentum() const
Returns 3-momentum at dca to primary vertex.
Definition: StMuTrack.h:260
StBTofNtupleMaker::Init
Int_t Init()
initialize ntuple and daqmap, and reset counters
Definition: StBTofNtupleMaker.cxx:94
kStOK
Definition: Stypes.h:40
StMuTrack::nHits
UShort_t nHits() const
Bingchu.
Definition: StMuTrack.h:237
StMuBTofPidTraits::timeOfFlight
float timeOfFlight() const
timing for PID
Definition: StMuBTofPidTraits.h:119
StMuTrack::dEdx
Double_t dEdx() const
Returns measured dE/dx value.
Definition: StMuTrack.h:248
StMuDst::btofHeader
static StBTofHeader * btofHeader()
returns pointer to the btofHeader - dongx
Definition: StMuDst.h:423
StMuTrack::nSigmaProton
Double_t nSigmaProton() const
Returns Craig&#39;s distance to the calculated dE/dx band for protons in units of sigma.
Definition: StMuTrack.h:247
StBTofNtupleMaker::Make
Int_t Make()
get tofr slat, pvpd rawdata and global data from StEvent and store in flat TTrees (ntuples) ...
Definition: StBTofNtupleMaker.cxx:195
StMuTrack::dcaGlobal
StThreeVectorF dcaGlobal(Int_t vtx_id=-1) const
Returns 3D distance of closest approach to primary vertex of associated global track.
Definition: StMuTrack.cxx:371
StMuTrack::nSigmaKaon
Double_t nSigmaKaon() const
Returns Craig&#39;s distance to the calculated dE/dx band for kaons in units of sigma.
Definition: StMuTrack.h:246
StBTofNtupleMaker::StBTofNtupleMaker
StBTofNtupleMaker(const Char_t *name, const Char_t *outname)
constructor sets default parameters
Definition: StBTofNtupleMaker.cxx:79
StBTofNtupleMaker::Finish
Int_t Finish()
write and close the ntuple file
Definition: StBTofNtupleMaker.cxx:172
Generated by   doxygen 1.8.5