StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StBTofCalibMaker.cxx
1 /*******************************************************************
2  *
3  * $Id: StBTofCalibMaker.cxx,v 1.24 2021/05/29 23:57:08 geurts Exp $
4  *
5  * Author: Xin Dong
6  *****************************************************************
7  *
8  * Description: - Tof Calibration Maker to do the calibration for pVPD
9  * (start timing) , TOF
10  * - store into StBTofPidTraits
11  *
12  *****************************************************************
13  *Revision 1.24 2021/12/06 15:11:00, bkimel
14  *select vertex in target region 198<V_z<202 when mFXTMode == true,
15  *include nT0 iterative outlier rejection for mFXTMode with new nT0 == 1 and nT0 == 2
16  *cases
17  *
18  *Revision 1.23 2020/10/09 11pm, Zaochen
19  *add (if (IAttr("btofFXT")) mFXTMode = kTRUE;) in the Init(),
20  *it could allow the chain option "btofFXT" to turn on the FXTMode easily
21 
22  *Revision 1.22 2020/04/09 4pm, Zaochen
23  *implement Xin's updates to allow more pions and protons for the T0s in FXT mode
24  *add a flag mFXTMode: 0 for Collider mode, 1 for FXT mode
25  *
26  * $Log: StBTofCalibMaker.cxx,v $
27  * Revision 1.24 2021/05/29 23:57:08 geurts
28  * Updates to improve pp and pA handling - by Bassam Aboona (TAMU)
29  *
30  * Revision 1.23 2021/01/27 04:06:25 geurts
31  * Introducing meaningful nTofSigma calculations in VPDstartless mode.
32  *
33  * Revision 1.22 2020/10/10 04:36:00 zye20
34  * new added chain option btofFXT which could turn on FXTMode of StBTofCalibMaker
35  *
36  * Revision 1.21 2020/10/10 04:31:03 zye20
37  * new added chain option btofFXT which could turn on FXTMode of StBTofCalibMaker
38  *
39  * Revision 1.20 2020/07/03 21:51:24 geurts
40  * Removed unnecessary warning messages and replaced them with counters.
41  *
42  * Revision 1.19 2020/04/10 20:41:38 zye20
43  * Xin add more pions and add protons for T0s in the FXT mode
44  *
45  * Revision 1.18 2019/04/23 05:49:57 jdb
46  * Added function to allow forcing 0 starttime for totally startless BTOF usage in UPC
47  *
48  * Revision 1.17 2017/10/20 17:50:32 smirnovd
49  * Squashed commit of the following:
50  *
51  * StBTof: Remove outdated ClassImp macro
52  *
53  * Prefer explicit namespace for std:: names in header files
54  *
55  * Removed unnecessary specification of default std::allocator
56  *
57  * Frank signed-off
58  *
59  * Revision 1.16 2017/03/02 18:30:44 jeromel
60  * Changes by jdb, nl - inData.open() of files on live disk TBF later
61  *
62  * Revision 1.10 2016/11/14 11:32:15 nluttrel
63  * Simulated hits no longer undergo electronics corrections
64  * If StVpdSimMaker used in chain, defaults to use Vpd start time
65  *
66  * Revision 1.15 2016/06/30 17:09:59 jdb
67  * Fixed Several errors identified by Coverity
68  *
69  * Revision 1.14 2011/07/27 15:44:32 geurts
70  * minor bug update: mProjVtxZ does not get initialized when mUseEventVertex is false, but is printed regardless.
71  *
72  * Revision 1.13 2011/05/09 14:32:10 geurts
73  * use appropriate log level for debug messages
74  *
75  * Revision 1.12 2010/10/31 05:52:11 geurts
76  * fixed module index range for read in loop for BOARD (TDIG) based calibration
77  *
78  * Revision 1.11 2010/10/30 05:20:50 geurts
79  * Calibration Maker reads (file/dbase) in and applies cell-based, module-based, or board-based (TDIG) calibration parameters
80  *
81  * Revision 1.10 2010/05/27 21:41:14 geurts
82  * Pick the default primary vertex (for mUseEventVertex). Additional cuts in selecting the vertex for tstart() have been removed.
83  *
84  * Revision 1.9 2010/05/25 22:09:18 geurts
85  * improved database handling and reduced log output
86  *
87  * Revision 1.8 2010/05/12 22:46:21 geurts
88  * Startless BTOF self-calibration method (Xin)
89  *
90  * Revision 1.7 2010/04/29 03:42:37 dongx
91  * Remove ranking>0 cut in event vertex selection for start time calculation
92  *
93  * Revision 1.6 2010/04/09 21:26:51 geurts
94  * Introduced "UseProjectedVertex" maker attribute to allow selection of the
95  * standard event vertex or one determined by track extrapolation
96  * (typically used in pp collisions).
97  *
98  * Revision 1.5 2010/04/03 15:43:58 dongx
99  * Change the default to use event vertex for start position for Run10 AuAu
100  *
101  * Revision 1.4 2010/03/04 23:10:20 dongx
102  * Added cleanup for PID variables in MuBTofPidTraits when processMuDst()
103  *
104  * Revision 1.3 2009/12/04 22:26:34 geurts
105  * Split original CalibMaker into dedicated StVpdCalibMaker and BTOF-specific StBTofCalibMaker (Xin):
106  * - function added to directly access the MuDst
107  * - clean up those VPD members and functions as they are moved to the StVpdCalibMaker
108  * - add VPD related functions to load/write the calibration VPD information in the BTofHeader
109  * - few small algorithm updates to be consistent with what is used in calibration procedures
110  * - several minor code cleanups
111  *
112  * Revision 1.2 2009/11/21 00:29:52 geurts
113  * Dtabase readout made more robust, static const moved to cxx.
114  *
115  * Revision 1.1 2009/09/23 02:28:41 geurts
116  * first version: Combined BTOF & VPD CalibMaker
117  *
118  *
119  *******************************************************************/
120 #include <iostream>
121 #include "StEvent/StEvent.h"
122 #include "StEvent/StBTofCollection.h"
123 #include "StEvent/StBTofHit.h"
124 #include "StEvent/StBTofHeader.h"
125 #include "StEvent/StBTofPidTraits.h"
126 #include "StEvent/StEventTypes.h"
127 #include "Stypes.h"
128 #include "StThreeVectorD.hh"
129 #include "StHelix.hh"
130 #include "StEvent/StTrackGeometry.h"
131 #include "StEvent/StTrackPidTraits.h"
132 #include "StEventUtilities/StuRefMult.hh"
133 #include "PhysicalConstants.h"
134 #include "StPhysicalHelixD.hh"
135 #include "tables/St_tofTOffset_Table.h"
136 #include "tables/St_tofTotbCorr_Table.h"
137 #include "tables/St_tofZbCorr_Table.h"
138 
139 #include "tables/St_vertexSeed_Table.h"
140 
141 #include "StBTofUtil/tofPathLength.hh"
142 #include "StBTofUtil/StBTofHitCollection.h"
143 #include "StBTofUtil/StBTofGeometry.h"
144 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
145 #include "StMuDSTMaker/COMMON/StMuDst.h"
146 #include "StMuDSTMaker/COMMON/StMuBTofHit.h"
147 #include "StMuDSTMaker/COMMON/StMuTrack.h"
148 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
149 #include "StMuDSTMaker/COMMON/StMuBTofPidTraits.h"
150 
151 #include "StBTofCalibMaker.h"
152 #include "StVpdCalibMaker/StVpdCalibMaker.h"
153 #include "StBTofUtil/StBTofSimResParams.h"
154 #include "StBTofUtil/StVpdSimConfig.h"
155 
156 #include "TProfile.h"
157 
159 const Double_t StBTofCalibMaker::VHRBIN2PS = 24.4140625; // 1000*25/1024 (ps/chn)
161 const Double_t StBTofCalibMaker::HRBIN2PS = 97.65625; // 97.65625= 1000*100/1024 (ps/chn)
163 const Double_t StBTofCalibMaker::TMAX = 51200.;
165 const Double_t StBTofCalibMaker::VZDIFFCUT=6.;
167 const Double_t StBTofCalibMaker::DCARCUT=1.;
168 const Double_t StBTofCalibMaker::mC_Light = C_C_LIGHT/1.e9;
169 
170 const Float_t StBTofCalibMaker::BTHBLCHCNST = 8.; // Bethe-Bloch constant used for dE/dx
171  // correction in start-time calculations
172  // in pppAMode
173 const Float_t StBTofCalibMaker::DEDXTCORR[2] = {0.033, 0.013}; // correction to the vpd hit times
174  // to match the corrected BTOF start time
175  // as a function of run year
176 
177 //_____________________________________________________________________________
178 StBTofCalibMaker::StBTofCalibMaker(const char *name) : StMaker(name)
179 {
183  setVPDHitsCut(1,1);
184  setOuterGeometry(true);
185 
186  mEvent = 0;
187  mBTofHeader = 0;
188  mMuDst = 0;
189  mZCalibType = NOTSET;
190  mTotCalibType = NOTSET;
191 
192  mSlewingCorr = kTRUE;
193  mMuDstIn = kFALSE;
194  mUseVpdStart = kTRUE;
195  mForceTStartZero = false;
196  isMcFlag = kFALSE;
197  mFXTMode = kFALSE;
198 
199  mPPPAMode = kFALSE;
200  mPPPAPionSel = kFALSE;
201  mPPPAOutlierRej = kFALSE;
202  mNSigmaTofMode = kFALSE;
203  mRun15Slew = kFALSE;
204  mPPPAModeHist = kFALSE;
205 
206  setCreateHistoFlag(kFALSE);
207  setHistoFileName("btofcalib.root");
208 
209  // default initialization from database
210  mInitFromFile = kFALSE;
211 
212  // assign default locations and names to the calibration files
213  setCalibFileTot("/star/institutions/rice/calib/default/totCali_4DB.dat");
214  setCalibFileZhit("/star/institutions/rice/calib/default/zCali_4DB.dat");
215  setCalibFileT0("/star/institutions/rice/calib/default/t0_4DB.dat");
216 
217 
218  StThreeVectorD MomFstPt(0.,0.,9999.);
219  StThreeVectorD origin(0.,0.,0.);
220  mBeamHelix = new StPhysicalHelixD(MomFstPt,origin,0.5*tesla,1.);
221 
222 
223 }
224 
225 //_____________________________________________________________________________
226 StBTofCalibMaker::~StBTofCalibMaker()
227 { /* noop */ }
228 
229 //_____________________________________________________________________________
230 void StBTofCalibMaker::resetPars()
231 {
232  memset(mTofTotEdge, 0, sizeof(mTofTotEdge));
233  memset(mTofTotCorr, 0, sizeof(mTofTotCorr));
234  memset(mTofZEdge, 0, sizeof(mTofZEdge) );
235  memset(mTofZCorr, 0, sizeof(mTofZCorr) );
236  memset(mTofTZero, 0, sizeof(mTofTZero) );
237 }
238 
239 //_____________________________________________________________________________
240 void StBTofCalibMaker::resetVpd()
241 {
242  memset(mVPDLeTime, 0, sizeof(mVPDLeTime));
243  mTStart = -9999.;
244  mTDiff = -9999.;
245  mProjVtxZ = -9999.;
246  mVPDVtxZ = -9999.;
247  mVPDHitPatternEast = 0;
248  mVPDHitPatternWest = 0;
249  mNEast = 0;
250  mNWest = 0;
251  mValidStartTime = kFALSE;
252  mNTzero = 0;
253 
254  mNTzeroCan = 0;
255  mTCanFirst = 99999.;
256  mTCanLast = -99999.;
257 
258  mVpdEHits = 0;
259  mVpdWHits = 0;
260  mVpdEGoodHits = 0;
261  mVpdWGoodHits = 0;
262  mEarliestVpdEHit = 99999.;
263  mEarliestVpdWHit = 99999.;
264  mClosestVpdEHit = 99999.;
265  mClosestVpdWHit = 99999.;
266  mLatestVpdEHit = -99999.;
267  mLatestVpdWHit = -99999.;
268 }
269 
270 //____________________________________________________________________________
271 Int_t StBTofCalibMaker::Init()
272 {
273  resetPars();
274  resetVpd();
275 
276 
277  if (IAttr("btofFXT")){//True for FXT mode calib, default as false for collider mode calib
278  mFXTMode = kTRUE;
279  LOG_INFO << "btofFXT mode is on." << endm;
280  }
281 
282 
283  if (IAttr("pppAMode")) {
284  mPPPAMode = kTRUE;
285  mRun15Slew = kTRUE;
286  LOG_INFO << "pppAMode is on." << endm;
287  }
288 
289  if (IAttr("setPPPAOutlierRej")) mPPPAOutlierRej = kTRUE;
290 
291  mUseEventVertex = ! IAttr("UseProjectedVertex");
292  if (mUseEventVertex) {
293  LOG_INFO << "Use event vertex position." << endm;
294  } else {
295  LOG_INFO << "Use projected vertex position." << endm;
296  }
297 
298  // m_Mode can be set by SetMode() method
299  if(m_Mode) {
300  // setHistoFileName("btofcalib.root");
301  } else {
302  setHistoFileName("");
303  }
304 
305  if (mHisto){
306  bookHistograms();
307  LOG_INFO << "Histograms are booked" << endm;
308  if (mHistoFileName!="") {
309  LOG_INFO << "Histograms will be stored in " << mHistoFileName.c_str() << endm;
310  }
311  }
312 
313  if (mPPPAModeHist){
314  bookPPPAHistograms();
315  LOG_INFO << "pppAMode Histograms are booked!" << endm;
316  }
317 
318  return kStOK;
319 }
320 
321 //____________________________________________________________________________
322 Int_t StBTofCalibMaker::InitRun(int runnumber)
323 {
324  // tof run configurations
325 
327  Int_t val = initParameters(runnumber);
328  if(val==kStOK) {
329  mValidCalibPar = kTRUE;
330  LOG_DEBUG << "Initialized valid calibration parameters." << endm;
331  } else {
332  mValidCalibPar = kFALSE;
333  LOG_WARN << "No valid calibration parameters! " << endm;
334  }
335 
337  StVpdCalibMaker *vpdCalib = (StVpdCalibMaker *)GetMaker("vpdCalib");
338 
340  mVpdResConfig = new StVpdSimConfig;
341  mVpdResConfig->loadVpdSimParams(); // do i really need this?
342  mVpdRes = mVpdResConfig->getParams();
343  mBTofRes = new StBTofSimResParams;
344  mBTofRes->loadParams(runnumber); //zaochen
345 
346 
347  if(vpdCalib) {
348  mUseVpdStart = vpdCalib->useVpdStart();
349 
350  if (mUseVpdStart) {LOG_INFO << "Found VPD Calibration Maker: use vpd for start timing" << endm;}
351  else {LOG_INFO << "Found VPD Calibration Maker: vpd **NOT** used for start timing" << endm;}
352  } else {
353  mUseVpdStart = kFALSE;
354  LOG_INFO << "NO VPD Calibration Maker found: vpd **NOT** used for start timing" << endm;
355  }
356 
357  // Apply fudge factor (FF) for dE/dx corrections (Runs 14 - 18)
358  if(runnumber>=14350003 && runnumber<=17315001) {
359  iYr = 0;
360  }
361  else {
362  iYr = 1;
363  }
364 
366  if(!mUseVpdStart && !mUseEventVertex) {
367  LOG_WARN << " Try to run calibration without VPD as the start time and DON'T use the event vertex! Wrong command! Exit!" << endm;
368  return kStOK;
369  }
370 
371  return kStOK;
372 
373 }
374 
375 
376 
377 //_____________________________________________________________________________
378 Int_t StBTofCalibMaker::initParameters(int runnumber)
379 {
382 
383  if (mInitFromFile){
384  LOG_INFO << "Initializing calibration parameters from files" << endm;
385  ifstream inData;
386 
388  LOG_INFO << " - ToT : " << mCalibFileTot << endm;
389  inData.open(mCalibFileTot.c_str());
390  unsigned int trayId, moduleId, cellId, boardId; // Coverity - jdb
391  int nbin;
392  int iCalibType; //9=cell, 8=module, 7=board -- Inset the number at the top of the .dat
393  inData >> iCalibType;
394 
395  // move enumeration to include file
396  //enum calibtype {board=960, module=3840, cell=23040} calibType;
397  mTotCalibType = calibtype(iCalibType);
398 
399  // switch(CalibType){ //selecting the calibration parameter format
400  switch(mTotCalibType) {
401  /********************************************************/
402  //cell based format
403  case CELLCALIB:
404  for(int i=0;i<mNTray;i++) {
405  for(int j=0;j<mNModule;j++) {
406  for(int l=0;l<mNCell;l++){
407  inData>>trayId>>moduleId>>cellId;
408  inData>>nbin;
409 
410  // coverity - jdb
411  // protect agains overflow
412  if ( trayId - 1 >= mNTray || moduleId - 1 >= mNModule || cellId - 1 >= mNCell ){
413  LOG_ERROR << "OUT OF BOUNDS, trayId = " << trayId << ", moduleId = " << moduleId << ", cellId = " << cellId << endl;
414  continue;
415  }
416  if ( nbin >= mNBinMax || nbin < 0){
417  LOG_ERROR << "OUT OF BOUNDS, # of TOT bins = " << nbin << endl;
418  continue;
419  }
420 
421  for( int k=0; k <= nbin; k++ ) {
422  inData >> mTofTotEdge[trayId-1][moduleId-1][cellId-1][k];
423  }
424 
425  for(int k=0;k <= nbin;k++) {
426  inData >> mTofTotCorr[trayId-1][moduleId-1][cellId-1][k];
427 
428  if(k%10==0&&Debug()) {
429  LOG_DEBUG << " ijlk= " << i << " " << j << " " << l << " " << k << " tot " << mTofTotEdge[trayId-1][moduleId-1][cellId-1][k] << " corr " << mTofTotCorr[trayId-1][moduleId-1][cellId-1][k] << endm;
430  }
431  }
432  }//cell
433  }//module
434  }//tray
435  break;
436  /********************************************************/
437  // Module based format
438  case MODULECALIB: //module based
439  for(int i=0;i<mNTray;i++) {
440  for(int j=0;j<mNModule;j++) {
441  inData>>trayId>>moduleId;
442  cellId = 1;
443 
444  inData>>nbin;
445 
446  // coverity - jdb
447  // protect agains overflow
448  if ( trayId - 1 >= mNTray || moduleId - 1 >= mNModule || cellId - 1 >= mNCell ){
449  LOG_ERROR << "OUT OF BOUNDS, trayId = " << trayId << ", moduleId = " << moduleId << ", cellId = " << cellId << endl;
450  continue;
451  }
452  if ( nbin >= mNBinMax || nbin < 0){
453  LOG_ERROR << "OUT OF BOUNDS, # of TOT bins = " << nbin << endl;
454  continue;
455  }
456 
457  for(int k=0;k<=nbin;k++){
458  inData>>mTofTotEdge[trayId-1][moduleId-1][cellId-1][k];
459 
460  // fill all cells
461  for(int l=0; l < mNCell; l++){
462  mTofTotEdge[trayId-1][moduleId-1][l][k] = mTofTotEdge[trayId-1][moduleId-1][cellId-1][k];
463  }
464  }
465  for(int k=0;k<=nbin;k++) {
466  inData>>mTofTotCorr[trayId-1][moduleId-1][cellId-1][k];
467  // fill all cells
468  for(int l=0;l<mNCell;l++){
469  mTofTotCorr[trayId-1][moduleId-1][l][k] = mTofTotCorr[trayId-1][moduleId-1][cellId-1][k];
470  }
471 
472  if(k%10==0&&Debug()) {
473  LOG_DEBUG << " ijlk= " << i << " " << j << " " << 0 << " " << k << " tot " << mTofTotEdge[trayId-1][moduleId-1][0][k] << " corr " << mTofTotCorr[trayId-1][moduleId-1][0][k] << endm;
474  }
475  }
476 
477  } //module
478  } //tray
479  break;
480 
481  /********************************************************/
482  // Board based format
483  case BOARDCALIB: //tdig based
484  for(int i=0;i<mNTray;i++) {
485  for(int j=0;j<mNTDIG;j++) {
486  inData>>trayId>>boardId;
487  // first in each board
488  moduleId = (boardId-1)*4;
489  cellId = 1;
490 
491  inData>>nbin;
492 
493  // coverity - jdb
494  // protect agains overflow
495  if ( trayId - 1 >= mNTray || moduleId >= mNModule || cellId - 1 >= mNCell || moduleId + 3 >= mNModule ){
496  LOG_ERROR << "OUT OF BOUNDS, trayId = " << trayId << ", moduleId = " << moduleId << ", cellId = " << cellId << endl;
497  continue;
498  }
499  if ( nbin >= mNBinMax || nbin < 0){
500  LOG_ERROR << "OUT OF BOUNDS, # of TOT bins = " << nbin << endl;
501  continue;
502  }
503 
504  for(int k=0;k<=nbin;k++){
505  inData>>mTofTotEdge[trayId-1][moduleId][cellId-1][k];
506 
507  for(int m=0;m<4;m++){
508  for(int l=0;l<mNCell;l++){
509  mTofTotEdge[trayId-1][moduleId+m][l][k] = mTofTotEdge[trayId-1][moduleId][cellId-1][k];
510  }//cell
511  }//modules per board
512  }
513 
514  for(int k=0;k<=nbin;k++) {
515 
516  inData>>mTofTotCorr[trayId-1][moduleId][cellId-1][k];
517 
518  for(int m=0;m<4;m++){
519  for(int l=0;l<mNCell;l++){
520  mTofTotCorr[trayId-1][moduleId+m][l][k] = mTofTotCorr[trayId-1][moduleId][cellId-1][k];
521  }//cell
522  }//modules per board
523 
524  if(k%10==0&&Debug()) {
525  LOG_DEBUG << " ijlk= " << i << " " << j*4 << " " << 0 << " " << k << " tot " << mTofTotEdge[trayId-1][moduleId][cellId][k] << " corr " << mTofTotCorr[trayId-1][moduleId][cellId][k] << endm;
526  }
527  } // k
528  }//board
529  }//tray
530  break;
531 
532  default:
533  LOG_WARN << "Please check the top of your TOT.dat file for the Calibration format. 9=cell,8=module,7=tdig. Your's is : " << mTotCalibType << endl;
534 
535  }//switch
536 
537  inData.close();
538 
539 
540  /******************************************************************************************************************************************************/
542  LOG_INFO << " - Zhit : " << mCalibFileZhit << endm;
543  inData.open(mCalibFileZhit.c_str());
544 
545  inData>>iCalibType;
546  mZCalibType = calibtype(iCalibType);
547 
548  switch(mZCalibType){ //selecting the calibration parameter format
549  case CELLCALIB: //cell based format
550  for(int i=0;i<mNTray;i++) {
551  for(int j=0;j<mNModule;j++) {
552  for(int l=0;l<mNCell;l++){
553  inData>>trayId>>moduleId>>cellId;
554  inData>>nbin;
555 
556  // coverity - jdb
557  // protect agains overflow
558  if ( trayId - 1 >= mNTray || moduleId - 1 >= mNModule || cellId - 1 >= mNCell ){
559  LOG_ERROR << "OUT OF BOUNDS, trayId = " << trayId << ", moduleId = " << moduleId << ", cellId = " << cellId << endl;
560  continue;
561  }
562  if ( nbin >= mNBinMax || nbin < 0){
563  LOG_ERROR << "OUT OF BOUNDS, # of TOT bins = " << nbin << endl;
564  continue;
565  }
566 
567  for(int k=0;k<=nbin;k++) inData>>mTofZEdge[trayId-1][moduleId-1][cellId-1][k];
568 
569  for(int k=0;k<=nbin;k++) {
570  inData>>mTofZCorr[trayId-1][moduleId-1][cellId-1][k];
571  if(k%10==0&&Debug()) {
572  LOG_DEBUG << " ijlk= " << i << " " << j << " " << l << " " << k << " zEdge " << mTofZEdge[trayId-1][moduleId-1][cellId-1][k] << " corr " << mTofZCorr[trayId-1][moduleId-1][cellId-1][k] << endm;
573  }
574  }
575 
576  }//cell
577  }//module
578  }//tray
579  break;
580 
581  case MODULECALIB: //module based
582  for(int i=0;i<mNTray;i++) {
583  for(int j=0;j<mNModule;j++) {
584  inData>>trayId>>moduleId;
585  cellId = 1;
586  inData>>nbin;
587 
588  // coverity - jdb
589  // protect agains overflow
590  if ( trayId - 1 >= mNTray || moduleId - 1 >= mNModule || cellId - 1 >= mNCell ){
591  LOG_ERROR << "OUT OF BOUNDS, trayId = " << trayId << ", moduleId = " << moduleId << ", cellId = " << cellId << endl;
592  continue;
593  }
594  if ( nbin >= mNBinMax || nbin < 0){
595  LOG_ERROR << "OUT OF BOUNDS, # of TOT bins = " << nbin << endl;
596  continue;
597  }
598 
599  for(int k=0;k<=nbin;k++){
600  inData>>mTofZEdge[trayId-1][moduleId-1][cellId-1][k];
601  for(int l=0;l<mNCell;l++){
602  mTofZEdge[trayId-1][moduleId-1][l][k]=mTofZEdge[trayId-1][moduleId-1][cellId-1][k];
603  }
604  }
605  for(int k=0;k<=nbin;k++) {
606  inData>>mTofZCorr[trayId-1][moduleId-1][cellId-1][k];
607  for(int l=0;l<mNCell;l++){
608  mTofZCorr[trayId-1][moduleId-1][l][k]=mTofZCorr[trayId-1][moduleId-1][cellId-1][k];
609  }
610  if(k%10==0&&Debug()) {
611  LOG_DEBUG << " ijlk= " << i << " " << j << " " << cellId-1 << " " << k << " zEdge " << mTofZEdge[trayId-1][moduleId-1][cellId][k] << " corr " << mTofZCorr[trayId-1][moduleId-1][cellId][k] << endm;
612  }
613  }
614  }//module
615  }//tray
616  break;
617 
618  case BOARDCALIB: //tdig based
619  for(int i=0;i<mNTray;i++) {
620  for(int j=0;j<mNTDIG;j++) {
621  inData>>trayId>>boardId;
622  // first in each board
623  moduleId = (boardId-1)*4;
624  cellId = 1;
625 
626  inData>>nbin;
627 
628  // coverity - jdb
629  // protect agains overflow
630  if ( trayId - 1 >= mNTray || moduleId >= mNModule || cellId - 1 >= mNCell || moduleId + 3 >= mNModule ){
631  LOG_ERROR << "OUT OF BOUNDS, trayId = " << trayId << ", moduleId = " << moduleId << ", cellId = " << cellId << endl;
632  continue;
633  }
634  if ( nbin >= mNBinMax || nbin < 0){
635  LOG_ERROR << "OUT OF BOUNDS, # of TOT bins = " << nbin << endl;
636  continue;
637  }
638 
639  for(int k=0;k<=nbin;k++){
640  inData>>mTofZEdge[trayId-1][moduleId][cellId-1][k];
641  for(int m=0;m<4;m++){
642 
643  for(int l=0;l<mNCell;l++){
644  mTofZEdge[trayId-1][moduleId+m][l][k] = mTofZEdge[trayId-1][moduleId][cellId-1][k];
645  }//cell
646  }//modules per board
647  } // k
648 
649  for(int k=0;k<=nbin;k++) {
650  inData>>mTofZCorr[trayId-1][moduleId][cellId-1][k];
651  for(int m=0;m<4;m++){
652  for(int l=0;l<mNCell;l++){
653  mTofZCorr[trayId-1][moduleId+m][l][k] = mTofZCorr[trayId-1][moduleId][cellId-1][k];
654  }//cell
655  }//modules per board
656 
657  if(k%10==0&&Debug()) {
658  LOG_DEBUG << " ijlk= " << i << " " << moduleId << " " << cellId-1 << " " << k << " tot " << mTofZEdge[trayId-1][moduleId][cellId][k] << " corr " << mTofZCorr[trayId-1][moduleId][cellId][k] << endm;
659  }
660  }// k bin
661  }//board
662  }//tray
663  break;
664 
665 
666  default:
667  LOG_WARN << "Please check the top of your Zhit.dat file for the Calibration format. 9=cell,8=module,7=tdig. Your's is : " << mZCalibType << endl;
668  }//switch
669  inData.close();
670 
671  /******************************************************************************************************************************************************/
673  LOG_INFO << " - T0 : " << mCalibFileT0 << endm;
674  inData.open(mCalibFileT0.c_str());
675  //int moduleId, cellId;
676  for(int i=0;i<mNTray;i++) {
677  for(int j=0;j<mNModule;j++) {
678  for(int k=0;k<mNCell;k++) {
679  inData>>trayId>>moduleId>>cellId;
680 
681  // coverity - jdb
682  // protect agains overflow
683  if ( trayId - 1 >= mNTray || moduleId - 1 >= mNModule || cellId - 1 >= mNCell ){
684  LOG_ERROR << "OUT OF BOUNDS, trayId = " << trayId << ", moduleId = " << moduleId << ", cellId = " << cellId << endl;
685  continue;
686  }
687 
688  inData>>mTofTZero[trayId-1][moduleId-1][cellId-1];
689  int index = i*mNModule*mNCell+j*mNCell+k;
690 
691  if(index%100==0&&Debug()) {
692  LOG_DEBUG << " ijk= " << i << " " << j << " " << k << " t0 " << mTofTZero[trayId-1][moduleId-1][cellId-1] << endm;
693  }
694  }
695  }
696  }
697  inData.close();
698 
699  // end load from file
700  } else {
701 
703  LOG_INFO << "Initializing calibration parameters from database" << endm;
704 
705  // read tofTotbCorr table
706  TDataSet *mDbDataSet = GetDataBase("Calibrations/tof/tofTotbCorr");
707  St_tofTotbCorr* tofTotCorr = static_cast<St_tofTotbCorr*>(mDbDataSet->Find("tofTotbCorr"));
708  if(!tofTotCorr) {
709  LOG_ERROR << "unable to get tof TotbCorr table parameters" << endm;
710  // assert(tofTotCorr);
711  return kStErr;
712  }
713  tofTotbCorr_st* totCorr = static_cast<tofTotbCorr_st*>(tofTotCorr->GetArray());
714  Int_t numRows = tofTotCorr->GetNRows();
715 
716  if(numRows!=mNTray*mNTDIG && numRows!=mNTray*mNModule*mNCell && numRows!=mNTray*mNModule) {
717  LOG_WARN << " Mis-matched number of rows in tofTotbCorr table! " << numRows
718  << " (exp:" << mNTray*mNTDIG << " or " << mNTray*mNModule*mNCell << " or "<< mNTray*mNModule << ")" << endm;
719  }
720 
721  LOG_INFO << " Number of rows read in: " << numRows << " for ToT correction" << endm;
722 
723  // convert to calibtype
724  mTotCalibType = calibtype(numRows);
725 
726  switch(mTotCalibType){
727  case CELLCALIB://cell
728  for (Int_t i=0;i<numRows;i++) {
729  short trayId = totCorr[i].trayId;
730  short moduleId = totCorr[i].moduleId;
731  short cellId = totCorr[i].cellId;
732  //short boardId = (moduleId-1)/4+1; // used for trays
733 
734  LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
735  for(Int_t j=0;j<mNBinMax;j++) {
736  if(trayId>0&&trayId<=mNTray&&moduleId>0&&moduleId<=mNModule&&cellId>0&&cellId<=mNCell){ // trays
737  mTofTotEdge[trayId-1][moduleId-1][cellId-1][j] = totCorr[i].tot[j];
738  mTofTotCorr[trayId-1][moduleId-1][cellId-1][j] = totCorr[i].corr[j];
739  if(Debug()&&j%10==0) { LOG_DEBUG << " j=" << j << " tot " << mTofTotEdge[trayId-1][moduleId-1][cellId-1][j] << " corr " << mTofTotCorr[trayId-1][moduleId-1][cellId-1][j] << endm; }
740  }
741  } // end j 0->mNBinMax
742  } // end i 0->numRows
743  break;
744 
745  case MODULECALIB://module
746  for (Int_t i=0;i<numRows;i++) {
747  short trayId = totCorr[i].trayId;
748  short moduleId = totCorr[i].moduleId;
749  short cellId = totCorr[i].cellId;
750  //short boardId = (moduleId-1)/4+1; // used for trays
751 
752  LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
753  for(Int_t j=0;j<mNBinMax;j++) {
754  if(trayId>0&&trayId<=mNTray&&moduleId>0&&moduleId<=mNModule){ // trays
755  for(Int_t k=0;k<mNCell;k++){
756  mTofTotEdge[trayId-1][moduleId-1][cellId-1+k][j] = totCorr[i].tot[j];
757  mTofTotCorr[trayId-1][moduleId-1][cellId-1+k][j] = totCorr[i].corr[j];
758  if(Debug()&&j%10==0) { LOG_DEBUG << " j=" << j << " tot " << mTofTotEdge[trayId-1][moduleId-1][cellId-1+k][j] << " corr " << mTofTotCorr[trayId-1][moduleId-1][cellId-1+k][j] << endm; }
759  }//duplicating entries into each cell
760  }
761  } // end j 0->mNBinMax
762  } // end i 0->numRows
763  break;
764 
765  case BOARDCALIB://board
766  for (Int_t i=0;i<numRows;i++) {
767  short trayId = totCorr[i].trayId;
768  short moduleId = totCorr[i].moduleId;
769  short cellId = totCorr[i].cellId;
770  short boardId = (moduleId-1)/4+1; // used for trays
771 
772  LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
773  for(Int_t j=0;j<mNBinMax;j++) {
774  if(trayId>0&&trayId<=mNTray&&boardId>0&&boardId<=mNTDIG){ // trays
775  for(Int_t k=0; k<4;k++){
776  for(Int_t l=0;l<mNCell;l++){
777  mTofTotEdge[trayId-1][moduleId-1+k][cellId-1+l][j] = totCorr[i].tot[j];
778  mTofTotCorr[trayId-1][moduleId-1+k][cellId-1+l][j] = totCorr[i].corr[j];
779  if(Debug()&&j%10==0) { LOG_DEBUG << " j=" << j << " tot " << mTofTotEdge[trayId-1][moduleId-1+k][cellId-1+l][j] << " corr " << mTofTotCorr[trayId-1][moduleId-1+k][cellId-1+l][j] << endm; }
780  }//duplicating into cells
781  }//duplication into modules
782  }
783  } // end j 0->mNBinMax
784  } // end i 0->numRows
785  break;
786 
787  default:
788  LOG_WARN << "Number of rows in tofTotbCorr table mis-matched. "<<endl;
789  }//end of switch
790 
791 
792  // read tofZbCorr table
793  mDbDataSet = GetDataBase("Calibrations/tof/tofZbCorr");
794  St_tofZbCorr* tofZCorr = static_cast<St_tofZbCorr*>(mDbDataSet->Find("tofZbCorr"));
795  if(!tofZCorr) {
796  LOG_ERROR << "unable to get tof ZbCorr table parameters" << endm;
797  // assert(tofZCorr);
798  return kStErr;
799  }
800  tofZbCorr_st* zCorr = static_cast<tofZbCorr_st*>(tofZCorr->GetArray());
801  numRows = tofZCorr->GetNRows();
802 
803  if(numRows!=mNTray*mNTDIG && numRows!=mNTray*mNModule*mNCell && numRows !=mNTray*mNModule) {
804  LOG_WARN << " Mis-matched number of rows in tofZbCorr table! " << numRows
805  << " (exp:" << mNTray*mNTDIG << " or " << mNTray*mNModule*mNCell << " or "<< mNTray*mNModule << ")" << endm;
806  }
807  LOG_INFO << " Number of rows read in: " << numRows << " for Z correction" << endm;
808 
809  // convert to calibtype
810  mZCalibType = calibtype(numRows);
811 
812  switch(mZCalibType){
813  case CELLCALIB://cell
814  for (Int_t i=0;i<numRows;i++) {
815  short trayId = totCorr[i].trayId;
816  short moduleId = totCorr[i].moduleId;
817  short cellId = totCorr[i].cellId;
818  //short boardId = (moduleId-1)/4+1; // used for trays
819 
820  LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
821  for(Int_t j=0;j<mNBinMax;j++) {
822  if(trayId>0&&trayId<=mNTray&&moduleId>0&&moduleId<=mNModule&&cellId>0&&cellId<=mNCell) { // trays
823  mTofZEdge[trayId-1][moduleId-1][cellId-1][j] = zCorr[i].z[j];
824  mTofZCorr[trayId-1][moduleId-1][cellId-1][j] = zCorr[i].corr[j];
825  if(Debug()&&j%10==0) { LOG_DEBUG << " j=" << j << " tot " << mTofZEdge[trayId-1][moduleId-1][cellId-1][j] << " corr " << mTofZCorr[trayId-1][moduleId-1][cellId-1][j] << endm; }
826  }
827  } // end j 0->mNBinMax
828  } // end i 0->numRows
829  break;
830 
831  case MODULECALIB://module
832  for (Int_t i=0;i<numRows;i++) {
833  short trayId = totCorr[i].trayId;
834  short moduleId = totCorr[i].moduleId;
835  short cellId = totCorr[i].cellId;
836  short boardId = (moduleId-1)/4+1; // used for trays
837 
838  LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
839  for(Int_t j=0;j<mNBinMax;j++) {
840  if(trayId>0&&trayId<=mNTray&&boardId>0&&moduleId<=mNModule) { // trays
841  for(Int_t k=0;k<mNCell;k++){
842  mTofZEdge[trayId-1][moduleId-1][cellId-1+k][j] = zCorr[i].z[j];
843  mTofZCorr[trayId-1][moduleId-1][cellId-1+k][j] = zCorr[i].corr[j];
844  if(Debug()&&j%10==0) { LOG_DEBUG << " j=" << j << " tot " << mTofZEdge[trayId-1][moduleId-1][cellId-1+k][j] << " corr " << mTofZCorr[trayId-1][moduleId-1][cellId-1+k][j] << endm; }
845  }//duplicating info to all cells
846  }
847  } // end j 0->mNBinMax
848  } // end i 0->numRows
849  break;
850 
851  case BOARDCALIB://board
852  for (Int_t i=0;i<numRows;i++) {
853  short trayId = totCorr[i].trayId;
854  short moduleId = totCorr[i].moduleId;
855  short cellId = totCorr[i].cellId;
856  short boardId = (moduleId-1)/4+1; // used for trays
857 
858  LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
859  for(Int_t j=0;j<mNBinMax;j++) {
860  if(trayId>0&&trayId<=mNTray&&boardId>0&&boardId<=mNTDIG) { // trays
861  for(Int_t k=0;k<4;k++){
862  for(Int_t l=0;l<mNCell;l++){
863  mTofZEdge[trayId-1][moduleId-1+k][cellId-1+l][j] = zCorr[i].z[j];
864  mTofZCorr[trayId-1][moduleId-1+k][cellId-1+l][j] = zCorr[i].corr[j];
865  if(Debug()&&j%10==0) { LOG_DEBUG << " j=" << j << " tot " << mTofZEdge[trayId-1][moduleId-1+k][cellId-1+l][j] << " corr " << mTofZCorr[trayId-1][moduleId-1+k][cellId-1+l][j] << endm; }
866  }//duplicating info to cell lvl
867  }//duplicating info to module lvl
868  }
869  } // end j 0->mNBinMax
870  } // end i 0->numRows
871  break;
872 
873  default:
874  LOG_WARN << "Number of rows in tofZbCorr table mis-matched. "<<endl;
875  }//switch
876  // read tofTOffset table
877  mDbDataSet = GetDataBase("Calibrations/tof/tofTOffset");
878  St_tofTOffset* tofTOffset = static_cast<St_tofTOffset*>(mDbDataSet->Find("tofTOffset"));
879  if(!tofTOffset) {
880  LOG_ERROR << "unable to get tof TOffset table parameters" << endm;
881  // assert(tofTOffset);
882  return kStErr;
883  }
884  tofTOffset_st* tZero = static_cast<tofTOffset_st*>(tofTOffset->GetArray());
885  numRows = tofTOffset->GetNRows();
886 
887  LOG_DEBUG << " Number of rows read in: " << numRows << " for TOffset correction" << endm;
888 
889  if(numRows!=mNTray) {
890  LOG_WARN << " Mis-matched number of rows in tofTOffset table! " << numRows
891  << " (exp:" << mNTray << ")" << endm;
892  // return kStErr;
893  }
894  for (Int_t i=0;i<numRows;i++) {
895  short trayId = tZero[i].trayId;
896  LOG_DEBUG << " tray " << trayId << endm;
897 
898  if(trayId>0&&trayId<=mNTray) {
899  for(int j=0;j<mNTOF;j++) {
900  mTofTZero[trayId-1][j/6][j%6] = tZero[i].T0[j];
901  if(Debug()&&j%10==0) { LOG_DEBUG << " j=" << j << " T0 " << mTofTZero[trayId-1][j/6][j%6] << endm; }
902  }
903  }
904  }
905  }
906 
907 
908  // ========== Set Beam Line =====================
909  double x0 = 0.;
910  double y0 = 0.;
911  double dxdz = 0.;
912  double dydz = 0.;
913 
914  // Get Current Beam Line Constraint from database
915  TDataSet* dbDataSet = this->GetDataBase("Calibrations/rhic/vertexSeed");
916 
917  if (dbDataSet) {
918  vertexSeed_st* vSeed = ((St_vertexSeed*) (dbDataSet->FindObject("vertexSeed")))->GetTable();
919 
920  x0 = vSeed->x0;
921  y0 = vSeed->y0;
922  dxdz = vSeed->dxdz;
923  dydz = vSeed->dydz;
924  }
925  else {
926  LOG_WARN << "No database for beamline (Calibrations/rhic/vertexSeed)" << endm;
927  }
928 
929  LOG_INFO << "BeamLine Constraint: " << endm;
930  LOG_INFO << "x(z) = " << x0 << " + " << dxdz << " * z" << endm;
931  LOG_INFO << "y(z) = " << y0 << " + " << dydz << " * z" << endm;
932 
933  //**********
934  //beam line not be calibrated yet
935  //x0 shift by 0.5
936  //x0 = 0.5;
937  //**********
938  StThreeVectorD origin(x0,y0,0.0);
939  double pt = 88889999;
940  double nxy=::sqrt(dxdz*dxdz + dydz*dydz);
941  if(nxy<1.e-5){ // beam line _MUST_ be tilted
942  LOG_WARN << "Beam line must be tilted!" << endm;
943  nxy=dxdz=1.e-5;
944  }
945  double p0=pt/nxy;
946  double px = p0*dxdz;
947  double py = p0*dydz;
948  double pz = p0; // approximation: nx,ny<<0
949  StThreeVectorD MomFstPt(px*GeV, py*GeV, pz*GeV);
950  if(mBeamHelix) delete mBeamHelix;
951  mBeamHelix = new StPhysicalHelixD(MomFstPt,origin,0.5*tesla,1.);
952 
953 
954  return kStOK;
955 }
956 
957 //____________________________________________________________________________
958 Int_t StBTofCalibMaker::FinishRun(int runnumber)
959 {
960  if(mBeamHelix) delete mBeamHelix;
961  mBeamHelix = 0;
962 
963  if (mBTofRes){delete mBTofRes; mBTofRes = 0;}
964  if (mVpdResConfig) {delete mVpdResConfig; mVpdResConfig = 0;}
965 
966  return kStOK;
967 }
968 
969 //_____________________________________________________________________________
970 Int_t StBTofCalibMaker::Finish()
971 {
972  if (mHistoFileName!="") writeHistograms();
973  if (mPPPAModeHist) writePPPAHistograms();
974  return kStOK;
975 }
976 
977 //_____________________________________________________________________________
978 Int_t StBTofCalibMaker::Make()
979 {
980  LOG_DEBUG << "StBTofCalibMaker::Maker: starting ..." << endm;
981  if(!mValidCalibPar) {
982  LOG_WARN << "No valid calibration parameters. Skip ... " << endm;
983  return kStOK;
984  }
985 
986  initEvent();
987  resetVpd();
988  if(mUseVpdStart) {
989  loadVpdData();
990  }
991 
992  if(!mMuDstIn) processStEvent();
993  else processMuDst();
994 
995  archiveVpdHitInfo();
996  writeStartTime();
997 
998  return kStOK;
999 }
1000 
1001 //_____________________________________________________________________________
1002 void StBTofCalibMaker::processStEvent()
1003 {
1004 
1005  // event selection // no primary vertex required
1006  if( !mEvent ) {LOG_WARN << "No StEvent" << endm; return;}
1007  if (!mEvent->btofCollection()) {LOG_WARN << "No BTOFCollection" << endm; return;}
1008  if (!mEvent->btofCollection()->hitsPresent()) {LOG_WARN << "No hits present" << endm; return;}
1009 
1010  StBTofCollection *theTof = mEvent->btofCollection();
1011  StSPtrVecBTofHit &tofHits = theTof->tofHits();
1012  Int_t nhits = tofHits.size();
1013  LOG_INFO << " Fired TOF cells + upVPD tubes : " << nhits << endm;
1014 
1015  if(mUseVpdStart) {
1016 
1017  mEvtVtxZ = -9999.;
1018  mProjVtxZ = -9999.;
1019  float dcaRmin = 9999.;
1020  // float rankHmax = -1.e9;
1021 
1022  if(mUseEventVertex) {
1023  // ///
1024  // /// select the vertex with highest positive rank within the VPDVtxZ cut range
1025  // ///
1026  // int nVtx = mEvent->numberOfPrimaryVertices();
1027  // for(int i=0;i<nVtx;i++) {
1028  // StPrimaryVertex *pVtx = mEvent->primaryVertex(i);
1029  // if(!pVtx) continue;
1030  // // if(pVtx->ranking()<0) continue; //! select positive ranking vertex
1031  // if(fabs(pVtx->position().z())>200.) continue; //! within 200 cm
1032  // if(fabs(pVtx->position().z()-mVPDVtxZ)>VZDIFFCUT) continue; //! VPDVtxZ cut
1033  // if(pVtx->ranking()<rankHmax) continue;
1034  // mEvtVtxZ = pVtx->position().z();
1035  // rankHmax = pVtx->ranking();
1036  // }
1038  // (Future version should allow for non-default vertex selections)
1039  StPrimaryVertex *pVtx = mEvent->primaryVertex();
1040 
1041  if (pVtx){
1042  mEvtVtxZ = pVtx->position().z();
1043  }
1044  else {
1045  LOG_WARN << "No (default) primary vertex information for this (st-) event" << endm;
1046  };
1047 
1048  tstart(mEvtVtxZ, &mTStart, &mTDiff);
1049 
1050  } else {
1054  for(int i=0;i<nhits;i++) {
1055  StBTofHit *aHit = dynamic_cast<StBTofHit*>(tofHits[i]);
1056  if(!aHit) continue;
1057  int trayId = aHit->tray();
1058  if(trayId>0&&trayId<=mNTray) {
1059  StGlobalTrack *gTrack = dynamic_cast<StGlobalTrack*>(aHit->associatedTrack());
1060  if(!gTrack) continue;
1061  StTrackGeometry *theTrackGeometry = gTrack->geometry();
1062 
1063  StThreeVectorD tofPos = theTrackGeometry->helix().at(theTrackGeometry->helix().pathLengths(*mBeamHelix).first);
1064  StThreeVectorD beamPos = mBeamHelix->at(theTrackGeometry->helix().pathLengths(*mBeamHelix).second);
1065  StThreeVectorD dcatof = tofPos - beamPos;
1066 
1067  LOG_DEBUG<<" tofPos(x,y,z) = "<<tofPos.x()<<","<<tofPos.y()<<","<<tofPos.z()<<endm;
1068  LOG_DEBUG<<"beamPos(x,y,z) = "<<beamPos.x()<<","<<beamPos.y()<<","<<beamPos.z()<<endm;
1069  LOG_DEBUG<<" dca (x,y,z) = "<<dcatof.x()<<","<<dcatof.y()<<","<<dcatof.z()<<endm;
1070  LOG_DEBUG<<" 2D dca = "<<sqrt(pow(dcatof.x(),2)+pow(dcatof.y(),2))<<endm;
1071  LOG_DEBUG<<" 2D signed dca = "<<theTrackGeometry->helix().geometricSignedDistance(beamPos.x(),beamPos.y())<<endm;
1072 
1074  if(fabs(tofPos.z()-mVPDVtxZ)>VZDIFFCUT) continue;
1075 
1076  if(dcaRmin>dcatof.perp()) {
1077  mProjVtxZ = tofPos.z();
1078  dcaRmin = dcatof.perp();
1079  }
1080  } // end if
1081  } // end loop tofhits
1082 
1083  if(dcaRmin>DCARCUT) mProjVtxZ = -9999.; // beam line contrain
1084  tstart(mProjVtxZ, &mTStart, &mTDiff);
1085 
1086  } // end if (mUseEventVertex)
1087 
1088  } else { // Don't use VPD as the start time
1089 
1090  StPrimaryVertex *pVtx = mEvent->primaryVertex();
1091  if(mFXTMode){
1092  int nPrimaryVertices = mEvent->numberOfPrimaryVertices();
1093  for(int iVtx=0; iVtx<nPrimaryVertices; iVtx++){
1094  pVtx = mEvent->primaryVertex(iVtx);
1095  if(pVtx->position().z() > 198. && pVtx->position().z() < 202.) break;
1096  else if(iVtx == nPrimaryVertices - 1){
1097  LOG_WARN << " No FXT primary vertex within target ... bye-bye" << endm;
1098  return;
1099  }
1100  }
1101  }
1102  if(!pVtx) {
1103  LOG_WARN << " No primary vertex ... bye-bye" << endm;
1104  return;
1105  }
1106  mEvtVtxZ = pVtx->position().z();
1107 
1108  tstart_NoVpd(theTof, pVtx, &mTStart);
1109 
1110  } // end if(mUseVpdStart)
1111 
1112  LOG_INFO << "primVz = " << mEvtVtxZ << " projVz = " << mProjVtxZ << " vpdVz = " << mVPDVtxZ << endm;
1113  LOG_INFO << "Tstart = " << mTStart << " Tdiff = " << mTDiff << " NTzero = " << mNTzero << endm;
1114  LOG_INFO << "NWest = " << mNWest << " NEast = " << mNEast << " TdcSum West = " << mTSumWest << " East = " << mTSumEast << endm;
1115 
1116 
1117  if(mTStart<-1000.) {
1118  LOG_INFO << "No valid start time for this event. Skip ..." << endm;
1119  mValidStartTime = kFALSE;
1120  return;
1121  } else {
1122  mValidStartTime = kTRUE;
1123  }
1124 
1125  //---------------------------------------
1126  // BTof calibration
1127  //---------------------------------------
1128 
1129  // keep a few counters:
1130  int nohitfound(0), ismchit(0),trayoutofbounds(0),notrack(0),
1131  gnopidtraits(0),gtraitisfalse(0), pidtoffalse(0),calibfailed(0),
1132  pgnopidtraits(0), ptraitisfalse(0),tofisset(0), num_doPID(0);
1133 
1134  for(int i=0;i<nhits;i++) {
1135 
1136  StBTofHit *aHit = dynamic_cast<StBTofHit*>(tofHits[i]);
1137  if(!aHit){
1138  LOG_DEBUG << "No hit found in the BTof Calibration." << endm;
1139  nohitfound++; continue;
1140  }
1141 
1142  isMcFlag = kFALSE; // Check to see if the hit is simulated.
1143  if ( aHit->qaTruth() == 1 ) {
1144  isMcFlag = kTRUE; ismchit++;
1145  LOG_DEBUG << "Simulated hit." << endm;
1146  }
1147 
1148  int trayId = aHit->tray();
1149  if(trayId<=0 || trayId>mNTray) {
1150  LOG_DEBUG << "Tray out of bounds." << endm;
1151  trayoutofbounds++; continue;
1152  }
1153 
1154  StGlobalTrack *gTrack = dynamic_cast<StGlobalTrack*>(aHit->associatedTrack());
1155  if(!gTrack) {
1156  LOG_DEBUG << " No associated Track with this hit." << endm;
1157  notrack++; continue;
1158  }
1159  else LOG_DEBUG << "Track found" << endm;
1160 
1161  const StPtrVecTrackPidTraits& theTofPidTraits = gTrack->pidTraits(kTofId);
1162  if(!theTofPidTraits.size()) {
1163  LOG_DEBUG << "Tof Pid Traits not populated." << endm;
1164  gnopidtraits++; continue;
1165  }
1166 
1167  StTrackPidTraits *theSelectedTrait = theTofPidTraits[theTofPidTraits.size()-1];
1168  if(!theSelectedTrait) {
1169  LOG_DEBUG << "The Selected Trait is false." << endm;
1170  gtraitisfalse++; continue;
1171  }
1172 
1173  StBTofPidTraits *pidTof = dynamic_cast<StBTofPidTraits *>(theSelectedTrait);
1174  if(!pidTof) {
1175  LOG_DEBUG << "PID Tof reports false" << endm;
1176  pidtoffalse++; continue;
1177  }
1178 
1179  double tot = aHit->tot(); // ns
1180  double tdc = aHit->leadingEdgeTime();
1181  double tof = tdc - mTStart;
1182  Double_t zhit = pidTof->zLocal();
1183 
1184  int moduleChan = (aHit->module()-1)*6 + (aHit->cell()-1);
1185  Double_t tofcorr = tof;
1186  if (!isMcFlag) {
1187  tofcorr = tofAllCorr(tof, tot, zhit, trayId, moduleChan);
1188  if(tofcorr<0.) {
1189  //LOG_WARN << " Calibration failed! ... " << endm;
1190  calibfailed++; continue;
1191  }
1192  }
1193 
1194  pidTof->setTimeOfFlight((Float_t)tofcorr);
1195 
1197  StPrimaryTrack *pTrack = dynamic_cast<StPrimaryTrack *>(gTrack->node()->track(primary));
1198  StBTofPidTraits *ppidTof = 0;
1199  if(pTrack) {
1200  const StPtrVecTrackPidTraits& pTofPidTraits = pTrack->pidTraits(kTofId);
1201  if(!pTofPidTraits.size()) {
1202  LOG_DEBUG << "pTofPidTraits is false." << endm;
1203  pgnopidtraits++; continue;
1204  }
1205 
1206  StTrackPidTraits *pSelectedTrait = pTofPidTraits[pTofPidTraits.size()-1];
1207  if(!pSelectedTrait) {
1208  LOG_DEBUG << "pSelectedTrait is false." << endm;
1209  ptraitisfalse++; continue;
1210  }
1211 
1212  ppidTof = dynamic_cast<StBTofPidTraits *>(pSelectedTrait);
1213 
1214  if(ppidTof && mUseEventVertex) {
1215  ppidTof->setTimeOfFlight((Float_t)tofcorr);
1216  tofisset++; LOG_DEBUG << "Time of Flight set." << endm;
1217  }
1218  }
1219 
1221  Double_t L = -9999.;
1222  Double_t ptot = -9999.;
1223  Bool_t doPID = kFALSE;
1224  if(mUseEventVertex) {
1225  if(!pTrack) {
1226  LOG_DEBUG << " The associated track is not a primary one. Skip PID calculation! " << endm;
1227  } else {
1228  StTrackGeometry *theTrackGeometry = pTrack->geometry();
1229  const StVertex *thisVertex = pTrack->vertex();
1230  if(!thisVertex) {
1231  LOG_DEBUG << " The associated track is not coming from any vertex. Skip PID calculation! " << endm;
1232  } else {
1233  StThreeVectorF primPos = thisVertex->position();
1234  L = tofPathLength(&primPos, &pidTof->position(), theTrackGeometry->helix().curvature());
1235  ptot = pTrack->geometry()->momentum().mag();
1236  doPID = kTRUE;
1237  LOG_DEBUG << "Pathlength and ptot set." << endm;
1238  }
1239  }
1240 
1241  } else {
1242 
1243  StTrackGeometry *theTrackGeometry = gTrack->geometry();
1244  StThreeVectorD tofPos = theTrackGeometry->helix().at(theTrackGeometry->helix().pathLengths(*mBeamHelix).first);
1245  StThreeVectorD dcatof = tofPos - mBeamHelix->at(theTrackGeometry->helix().pathLengths(*mBeamHelix).second);
1246  if(dcatof.perp()>DCARCUT) {
1247  LOG_DEBUG << " The projected position is far from beam line. Skip PID calculation! " << endm;
1248  } else if(fabs(tofPos.z()-mVPDVtxZ)>VZDIFFCUT) {
1249  LOG_DEBUG << " This track is not coming from the same VPD vertex! Skip PID calculation! " << endm;
1250  } else {
1251  L = tofPathLength(&tofPos, &pidTof->position(), theTrackGeometry->helix().curvature());
1252  ptot = gTrack->geometry()->momentum().mag();
1253  if(gTrack->dcaGeometry()) {
1254  ptot = gTrack->dcaGeometry()->momentum().mag();
1255  }
1256  doPID = kTRUE;
1257  }
1258 
1259  }
1260 
1261  if(!doPID) continue;
1262  num_doPID++;
1263 
1264  Double_t beta = L/(tofcorr*(C_C_LIGHT/1.e9));
1265 
1266  Double_t b_e = ptot/sqrt(ptot*ptot+M_ELECTRON*M_ELECTRON);
1267  Double_t b_pi = ptot/sqrt(ptot*ptot+M_PION_PLUS*M_PION_PLUS);
1268  Double_t b_k = ptot/sqrt(ptot*ptot+M_KAON_PLUS*M_KAON_PLUS);
1269  Double_t b_p = ptot/sqrt(ptot*ptot+M_PROTON*M_PROTON);
1270 
1271  float sigmae = -9999.;
1272  float sigmapi = -9999.;
1273  float sigmak = -9999.;
1274  float sigmap = -9999.;
1275 // float res = 0.013; // 0.013 by default - 1/beta resolution
1276  float res = tofCellResolution(trayId, moduleChan);
1277  float res_c = res * (C_C_LIGHT/1.e9);
1278 
1279  if(fabs(res)>1.e-5) {
1280  sigmae = (Float_t)(L*(1./beta-1./b_e)/res_c);
1281  sigmapi = (Float_t)(L*(1./beta-1./b_pi)/res_c);
1282  sigmak = (Float_t)(L*(1./beta-1./b_k)/res_c);
1283  sigmap = (Float_t)(L*(1./beta-1./b_p)/res_c);
1284  }
1285 
1286  pidTof->setPathLength((Float_t)L);
1287  pidTof->setBeta((Float_t)beta);
1288  pidTof->setSigmaElectron(sigmae);
1289  pidTof->setSigmaPion(sigmapi);
1290  pidTof->setSigmaKaon(sigmak);
1291  pidTof->setSigmaProton(sigmap);
1292 
1293  LOG_DEBUG << " storing BTofPidTraits for the global track" << endm;
1294 
1295  if(mUseEventVertex) {
1296 
1297  if(ppidTof) {
1298 
1299  ppidTof->setPathLength((Float_t)L);
1300  ppidTof->setBeta((Float_t)beta);
1301  ppidTof->setSigmaElectron(sigmae);
1302  ppidTof->setSigmaPion(sigmapi);
1303  ppidTof->setSigmaKaon(sigmak);
1304  ppidTof->setSigmaProton(sigmap);
1305 
1306  LOG_DEBUG << " storing BTofPidTraits for the primary track" << endm;
1307  } // end if ppidTof
1308  } // end if mUseEventVertex
1309 
1310  } // end tof hits
1311 
1312  LOG_INFO << "nohitfound:"<< nohitfound << " ismchit:" <<ismchit << " trayoutofbounds:" << trayoutofbounds << " notrack:" << notrack << endm;
1313  LOG_INFO << " gnopidtraits:" << gnopidtraits <<" gtraitisfalse:" << gtraitisfalse << " pidtoffalse:"<< pidtoffalse <<" calibfailed:"<< calibfailed << endm;
1314  LOG_INFO << " pgnopidtraits:"<< pgnopidtraits << " ptraitisfalse:" << ptraitisfalse << " tofisset:"<< tofisset << " doPID:" << num_doPID << endm;
1315 
1316  return;
1317 }
1318 
1319 //_____________________________________________________________________________
1320 void StBTofCalibMaker::processMuDst()
1321 {
1322  if(!mMuDst) {
1323  LOG_WARN << " No MuDst ... bye-bye" << endm;
1324  return;
1325  }
1326 
1327  cleanCalibMuDst();
1328 
1329  Int_t nhits = mMuDst->numberOfBTofHit();
1330  LOG_INFO << " Fired TOF cells + upVPD tubes : " << nhits << endm;
1331 
1332  if(mUseVpdStart) {
1333 
1334  mEvtVtxZ = -9999.;
1335  mProjVtxZ = -9999.;
1336  float dcaRmin = 9999.;
1337  // float rankHmax = -1.e9;
1338 
1339  if(mUseEventVertex) {
1340  // ///
1341  // /// select the vertex with highest positive rank within the VPDVtxZ cut range
1342  // ///
1343  // int nVtx = mMuDst->numberOfPrimaryVertices();
1344  // for(int i=0;i<nVtx;i++) {
1345  // StMuPrimaryVertex* pVtx = mMuDst->primaryVertex(i);
1346  // if(!pVtx) continue;
1347  // // if(pVtx->ranking()<0) continue; //! select positive ranking vertex
1348  // if(fabs(pVtx->position().z())>200.) continue; //! within 200 cm
1349  // if(fabs(pVtx->position().z()-mVPDVtxZ)>VZDIFFCUT) continue; //! VPDVtxZ cut
1350  // if(pVtx->ranking()<rankHmax) continue;
1351  // mEvtVtxZ = pVtx->position().z();
1352  // rankHmax = pVtx->ranking();
1353  // }
1355  // (Future version should allow for non-default vertex selections)
1356  StMuPrimaryVertex* pVtx = mMuDst->primaryVertex();
1357  if (pVtx){
1358  mEvtVtxZ = pVtx->position().z();
1359  }
1360  else {
1361  LOG_WARN << "No (default) primary vertex information for this (mudst) event" << endm;
1362  }
1363 
1364  tstart(mEvtVtxZ, &mTStart, &mTDiff);
1365 
1366  } else {
1370  for(int i=0;i<nhits;i++) {
1371  StMuBTofHit *aHit = (StMuBTofHit*)mMuDst->btofHit(i);
1372  if(!aHit) continue;
1373  int trayId = aHit->tray();
1374  if(trayId>0&&trayId<=mNTray) {
1375  StMuTrack *gTrack = aHit->globalTrack();
1376  if(!gTrack) continue;
1377 
1378  StPhysicalHelixD thisHelix = gTrack->helix();
1379 
1380  StThreeVectorD tofPos = thisHelix.at(thisHelix.pathLengths(*mBeamHelix).first);
1381  StThreeVectorD dcatof = tofPos - mBeamHelix->at(thisHelix.pathLengths(*mBeamHelix).second);
1382 
1384  if(fabs(tofPos.z()-mVPDVtxZ)>VZDIFFCUT) continue;
1385 
1386  if(dcaRmin>dcatof.perp()) {
1387  mProjVtxZ = tofPos.z();
1388  dcaRmin = dcatof.perp();
1389  }
1390  } // end if
1391  } // end loop tofhits
1392 
1393  if(dcaRmin>DCARCUT) mProjVtxZ = -9999.; // beam line contrain
1394  tstart(mProjVtxZ, &mTStart, &mTDiff);
1395 
1396  } // end if(mUseEventVertex)
1397  } else { // don't use vpd as the start time
1398 
1399  StMuPrimaryVertex *pVtx = mMuDst->primaryVertex();
1400  if(mFXTMode){
1401  int nPrimaryVertices = mMuDst->numberOfPrimaryVertices();
1402  for(int iVtx=0; iVtx<nPrimaryVertices; iVtx++){
1403  pVtx = mMuDst->primaryVertex(iVtx);
1404  if(pVtx->position().z() > 198. && pVtx->position().z() < 202.) break;
1405  else if(iVtx == nPrimaryVertices - 1){
1406  LOG_WARN << " No FXT primary vertex within target ... bye-bye" <<endm;
1407  return;
1408  }
1409  }
1410  }
1411  if(!pVtx) {
1412  LOG_WARN << " No primary vertex ... bye-bye" << endm;
1413  return;
1414  }
1415  mEvtVtxZ = pVtx->position().z();
1416 
1417  tstart_NoVpd(mMuDst, pVtx, &mTStart);
1418  }
1419 
1420  LOG_INFO << "primVz = " << mEvtVtxZ << " projVz = " << mProjVtxZ << " vpdVz = " << mVPDVtxZ << endm;
1421  LOG_INFO << "Tstart = " << mTStart << " Tdiff = " << mTDiff << " NTzero = " << mNTzero << endm;
1422  LOG_INFO << "NWest = " << mNWest << " NEast = " << mNEast << " TdcSum West = " << mTSumWest << " East = " << mTSumEast << endm;
1423 
1424 
1425  if(mTStart<-1000.) {
1426  LOG_INFO << " No valid start time for this event. Skip ..." << endm;
1427  mValidStartTime = kFALSE;
1428  return;
1429  } else {
1430  mValidStartTime = kTRUE;
1431  }
1432 
1433  //---------------------------------------
1434  // BTof calibration
1435  //---------------------------------------
1436 
1437  // keep a few counters:
1438  int nohitfound(0), ismchit(0),trayoutofbounds(0),notrack(0),
1439  gnopidtraits(0),gtraitisfalse(0), pidtoffalse(0),calibfailed(0),
1440  pgnopidtraits(0), ptraitisfalse(0),tofisset(0), num_doPID(0);
1441 
1442  for(int i=0;i<nhits;i++) {
1443  StMuBTofHit *aHit = (StMuBTofHit*)mMuDst->btofHit(i);
1444  if(!aHit) {nohitfound++; continue;}
1445  int trayId = aHit->tray();
1446  if(trayId<=0 || trayId>mNTray) {trayoutofbounds++;continue;}
1447 
1448  StMuTrack *gTrack = aHit->globalTrack();
1449  if(!gTrack) {
1450  LOG_DEBUG << " No associated Track with this hit." << endm;
1451  notrack++; continue;
1452  }
1453 
1454  isMcFlag = kFALSE; // Check to see if the hit is simulated.
1455  if ( aHit->qaTruth() == 1 ) {
1456  isMcFlag = kTRUE; ismchit++;
1457  LOG_DEBUG << "Simulated hit." << endm;
1458  }
1459 
1460  StMuBTofPidTraits pidTof = gTrack->btofPidTraits();
1461 
1462  double tot = aHit->tot(); // ns
1463  double tdc = aHit->leadingEdgeTime();
1464  while(tdc>TMAX) tdc -= TMAX;
1465  double tof = tdc - mTStart;
1466  Double_t zhit = pidTof.zLocal();
1467 
1468  int moduleChan = (aHit->module()-1)*6 + (aHit->cell()-1);
1469 
1470  Double_t tofcorr = tof;
1471  if (!isMcFlag) {
1472  tofcorr = tofAllCorr(tof, tot, zhit, trayId, moduleChan);
1473  if(tofcorr<0.) {
1474  //LOG_WARN << " Calibration failed! ... " << endm;
1475  calibfailed++; continue;
1476  }
1477  }
1478 
1480  pidTof.setTimeOfFlight((Float_t)tofcorr);
1481 
1483  StMuTrack *pTrack = aHit->primaryTrack();
1484  StMuBTofPidTraits ppidTof;
1485  if(pTrack) {
1486  ppidTof = pTrack->btofPidTraits();
1487  if(mUseEventVertex) ppidTof.setTimeOfFlight((Float_t)tofcorr);
1488  }
1489 
1491  Double_t L = -9999.;
1492  Double_t ptot = -9999.;
1493  Bool_t doPID = kFALSE;
1494  if(mUseEventVertex) {
1495  if(!pTrack) {
1496  LOG_DEBUG << " The associated track is not a primary one. Skip PID calculation! " << endm;
1497  } else {
1498  int iv = pTrack->vertexIndex();
1499  StMuPrimaryVertex *thisVertex = mMuDst->primaryVertex(iv);
1500  if(!thisVertex) {
1501  LOG_DEBUG << " The associated track is not coming from any vertex. Skip PID calculation! " << endm;
1502  } else {
1503  StThreeVectorF primPos = thisVertex->position();
1504  StPhysicalHelixD thisHelix = pTrack->helix();
1505  L = tofPathLength(&primPos, &pidTof.position(), thisHelix.curvature());
1506  ptot = pTrack->momentum().mag();
1507  doPID = kTRUE;
1508  }
1509  }
1510 
1511  } else {
1512 
1513  StPhysicalHelixD gHelix = gTrack->helix();
1514  StThreeVectorD tofPos = gHelix.at(gHelix.pathLengths(*mBeamHelix).first);
1515  StThreeVectorD dcatof = tofPos - mBeamHelix->at(gHelix.pathLengths(*mBeamHelix).second);
1516  if(dcatof.perp()>DCARCUT) {
1517  LOG_DEBUG << " The projected position is far from beam line. Skip PID calculation! " << endm;
1518  } else if(fabs(tofPos.z()-mVPDVtxZ)>VZDIFFCUT) {
1519  LOG_DEBUG << " This track is not coming from the same VPD vertex! Skip PID calculation! " << endm;
1520  } else {
1521  L = tofPathLength(&tofPos, &pidTof.position(), gHelix.curvature());
1522  ptot = gTrack->momentum().mag();
1523  doPID = kTRUE;
1524  }
1525  }
1526 
1527  if(doPID) {
1528  num_doPID++;
1529  Double_t beta = L/(tofcorr*(C_C_LIGHT/1.e9));
1530 
1531  Double_t b_e = ptot/sqrt(ptot*ptot+M_ELECTRON*M_ELECTRON);
1532  Double_t b_pi = ptot/sqrt(ptot*ptot+M_PION_PLUS*M_PION_PLUS);
1533  Double_t b_k = ptot/sqrt(ptot*ptot+M_KAON_PLUS*M_KAON_PLUS);
1534  Double_t b_p = ptot/sqrt(ptot*ptot+M_PROTON*M_PROTON);
1535 
1536  float sigmae = -9999.;
1537  float sigmapi = -9999.;
1538  float sigmak = -9999.;
1539  float sigmap = -9999.;
1540 // float res = 0.013; // 0.013 by default - 1/beta resolution
1541  float res = tofCellResolution(trayId, moduleChan);
1542  float res_c = res * (C_C_LIGHT/1.e9);
1543 
1544  if(fabs(res)>1.e-5) {
1545  sigmae = (Float_t)(L*(1./beta-1./b_e)/res_c);
1546  sigmapi = (Float_t)(L*(1./beta-1./b_pi)/res_c);
1547  sigmak = (Float_t)(L*(1./beta-1./b_k)/res_c);
1548  sigmap = (Float_t)(L*(1./beta-1./b_p)/res_c);
1549  }
1550 
1551  pidTof.setPathLength((Float_t)L);
1552  pidTof.setBeta((Float_t)beta);
1553  pidTof.setSigmaElectron(sigmae);
1554  pidTof.setSigmaPion(sigmapi);
1555  pidTof.setSigmaKaon(sigmak);
1556  pidTof.setSigmaProton(sigmap);
1557 
1558  if(mUseEventVertex && pTrack) {
1559 
1560  ppidTof.setPathLength((Float_t)L);
1561  ppidTof.setBeta((Float_t)beta);
1562  ppidTof.setSigmaElectron(sigmae);
1563  ppidTof.setSigmaPion(sigmapi);
1564  ppidTof.setSigmaKaon(sigmak);
1565  ppidTof.setSigmaProton(sigmap);
1566  }
1567  }
1568 
1569  gTrack->setBTofPidTraits(pidTof);
1570  LOG_DEBUG << " storing BTofPidTraits for the global track" << endm;
1571 
1572  if(mUseEventVertex && pTrack) {
1573  pTrack->setBTofPidTraits(ppidTof);
1574  LOG_DEBUG << " storing BTofPidTraits for the primary track" << endm;
1575  }
1576  } // end tof hits
1577 
1578  LOG_INFO << "nohitfound:"<< nohitfound << " ismchit:" <<ismchit << " trayoutofbounds:" << trayoutofbounds << " notrack:" << notrack << endm;
1579  LOG_INFO << " gnopidtraits:" << gnopidtraits <<" gtraitisfalse:" << gtraitisfalse << " pidtoffalse:"<< pidtoffalse <<" calibfailed:"<< calibfailed << endm;
1580  LOG_INFO << " pgnopidtraits:"<< pgnopidtraits << " ptraitisfalse:" << ptraitisfalse << " tofisset:"<< tofisset << " doPID:" << num_doPID << endm;
1581 
1582  return;
1583 }
1584 
1585 
1586 //_____________________________________________________________________________
1587 void StBTofCalibMaker::cleanCalibMuDst()
1588 {
1589  if(!mMuDst) return;
1590 
1591  Int_t nPrimary = mMuDst->numberOfPrimaryTracks();
1592  Int_t nGlobal = mMuDst->numberOfGlobalTracks();
1593  for(int i=0;i<nPrimary;i++) {
1594  StMuTrack *pTrack = (StMuTrack *)mMuDst->primaryTracks(i);
1595  if(!pTrack) continue;
1596  StMuBTofPidTraits pid = pTrack->btofPidTraits();
1597  cleanCalib(pid);
1598  pTrack->setBTofPidTraits(pid);
1599  }
1600  for(int i=0;i<nGlobal;i++) {
1601  StMuTrack *gTrack = (StMuTrack *)mMuDst->globalTracks(i);
1602  if(!gTrack) continue;
1603  StMuBTofPidTraits pid = gTrack->btofPidTraits();
1604  cleanCalib(pid);
1605  gTrack->setBTofPidTraits(pid);
1606  }
1607  return;
1608 }
1609 
1610 //_____________________________________________________________________________
1611 void StBTofCalibMaker::cleanCalib(StMuBTofPidTraits& pid)
1612 {
1613  pid.setTimeOfFlight(-999.);
1614  pid.setPathLength(-999.);
1615  pid.setBeta(-999.);
1616  pid.setSigmaElectron(-999.);
1617  pid.setSigmaPion(-999.);
1618  pid.setSigmaKaon(-999.);
1619  pid.setSigmaProton(-999.);
1620  pid.setProbElectron(-999.);
1621  pid.setProbPion(-999.);
1622  pid.setProbKaon(-999.);
1623  pid.setProbProton(-999.);
1624  return;
1625 }
1626 
1627 //_____________________________________________________________________________
1628 void StBTofCalibMaker::initEvent()
1629 {
1630  if(mMuDstIn) {
1631  StMuDstMaker *mMuDstMaker = (StMuDstMaker *)GetMaker("MuDst");
1632  if(!mMuDstMaker) {
1633  LOG_WARN << " No MuDstMaker ... bye-bye" << endm;
1634  return;
1635  }
1636  mMuDst = mMuDstMaker->muDst();
1637  if(!mMuDst) {
1638  LOG_WARN << " No MuDst ... bye-bye" << endm;
1639  return;
1640  }
1641 
1642  mBTofHeader = mMuDst->btofHeader();
1643  } else {
1644  mEvent = (StEvent *) GetInputDS("StEvent");
1645 
1646  if(!mEvent) return;
1647  StBTofCollection *btofColl = mEvent->btofCollection();
1648  if(!btofColl) return;
1649  mBTofHeader = btofColl->tofHeader();
1650  }
1651 
1652  return;
1653 }
1654 
1655 //_____________________________________________________________________________
1656 void StBTofCalibMaker::loadVpdData()
1657 {
1658  if(!mBTofHeader) return;
1659 
1660  mTSumWest = 0;
1661  mTSumEast = 0;
1662  mTSumWestSigma = 0;
1663  mTSumEastSigma = 0;
1664  mVPDHitPatternWest = mBTofHeader->vpdHitPattern(west);
1665  mVPDHitPatternEast = mBTofHeader->vpdHitPattern(east);
1666  mNWest = mBTofHeader->numberOfVpdHits(west);
1667  mNEast = mBTofHeader->numberOfVpdHits(east);
1668  mVPDVtxZ = mBTofHeader->vpdVz();
1669 
1670  for(int i=0;i<mNVPD;i++) {
1671  mVPDLeTime[i] = mBTofHeader->vpdTime(west, i+1);
1672  if(mVPDLeTime[i]>0.) {
1673  mTSumWest += mVPDLeTime[i];
1674  //fg add here the mTSumEastSigma based on what tubes were used.
1675  }
1676  if(Debug()) {
1677  LOG_DEBUG << " loading VPD West tubeId = " << i+1 << " time = " << mVPDLeTime[i] << endm;
1678  }
1679  }
1680 
1681  for(int i=0;i<mNVPD;i++) {
1682  mVPDLeTime[i+mNVPD] = mBTofHeader->vpdTime(east, i+1);
1683  if(mVPDLeTime[i+mNVPD]>0.) mTSumEast += mVPDLeTime[i+mNVPD];
1684  if(Debug()) {
1685  LOG_DEBUG << " loading VPD East tubeId = " << i+1 << " time = " << mVPDLeTime[i+mNVPD] << endm;
1686  }
1687  }
1688 
1689  return;
1690 }
1691 
1692 //_____________________________________________________________________________
1693 void StBTofCalibMaker::writeStartTime()
1694 {
1695  if(mBTofHeader) {
1696  mBTofHeader->setTStart(mTStart);
1697  mBTofHeader->setTDiff(mTDiff);
1698  mBTofHeader->setNTzero(mNTzero);
1699 
1700  mBTofHeader->setNTzeroCan(mNTzeroCan);
1701  mBTofHeader->setTCanFirst(mTCanFirst);
1702  mBTofHeader->setTCanLast(mTCanLast);
1703 
1704  mBTofHeader->setVpdEHits(mVpdEHits);
1705  mBTofHeader->setVpdWHits(mVpdWHits);
1706  mBTofHeader->setVpdEGoodHits(mVpdEGoodHits);
1707  mBTofHeader->setVpdWGoodHits(mVpdWGoodHits);
1708  mBTofHeader->setEarliestVpdEHit(mEarliestVpdEHit);
1709  mBTofHeader->setEarliestVpdWHit(mEarliestVpdWHit);
1710  mBTofHeader->setClosestVpdEHit(mClosestVpdEHit);
1711  mBTofHeader->setClosestVpdWHit(mClosestVpdWHit);
1712  mBTofHeader->setLatestVpdEHit(mLatestVpdEHit);
1713  mBTofHeader->setLatestVpdWHit(mLatestVpdWHit);
1714  }
1715 
1716  return;
1717 }
1718 
1719 //_____________________________________________________________________________
1720 Double_t StBTofCalibMaker::tofAllCorr(const Double_t tof, const Double_t tot, const Double_t z, const Int_t iTray, const Int_t iModuleChan)
1721 {
1722  int tray = iTray;
1723  int module = iModuleChan/6 + 1;
1724  int cell = iModuleChan%6 + 1;
1725  //int board = iModuleChan/24 + 1;
1726  LOG_DEBUG << "\nStBTofCalibMaker::btofAllCorr: BTof calibrating...\n"
1727  << "\tDoing Calibration in BTOF Tray " << tray << " Module " << module << " Cell " << cell
1728  << "\n\tinput tof = " << tof
1729  << " TOT = " << tot << " Zlocal = " << z << endm;
1730 
1731  Double_t tofcorr = tof;
1732 
1733  tofcorr -= mTofTZero[tray-1][module-1][cell-1];
1734 
1735  LOG_DEBUG << "T0 correction: "<<mTofTZero[tray-1][module-1][cell-1]<<endm;
1736 
1737  if(mSlewingCorr) {
1738  int iTotBin = -1;
1739  for(int i=0;i<mNBinMax-1;i++) {
1740  if(tot>=mTofTotEdge[tray-1][module-1][cell-1][i] && tot<mTofTotEdge[tray-1][module-1][cell-1][i+1]) {
1741  iTotBin = i;
1742  break;
1743  }
1744  }
1745  if(iTotBin>=0&&iTotBin<mNBinMax) {
1746  double x1 = mTofTotEdge[tray-1][module-1][cell-1][iTotBin];
1747  double x2 = mTofTotEdge[tray-1][module-1][cell-1][iTotBin+1];
1748  double y1 = mTofTotCorr[tray-1][module-1][cell-1][iTotBin];
1749  double y2 = mTofTotCorr[tray-1][module-1][cell-1][iTotBin+1];
1750  double dcorr = y1 + (tot-x1)*(y2-y1)/(x2-x1);
1751  LOG_DEBUG << "TOT correction: "<<dcorr<<endm;
1752 
1753  tofcorr -= dcorr;
1754  } else {
1755  LOG_DEBUG << " TOT out of range! EXIT! " << endm;
1756  return -9999.;
1757  }
1758 
1759  int iZBin = -1;
1760  for(int i=0;i<mNBinMax-1;i++) {
1761  if(z>=mTofZEdge[tray-1][module-1][cell-1][i] && z<mTofZEdge[tray-1][module-1][cell-1][i+1]) {
1762  iZBin = i;
1763  break;
1764  }
1765  }
1766  if(iZBin>=0&&iZBin<mNBinMax) {
1767  double x1 = mTofZEdge[tray-1][module-1][cell-1][iZBin];
1768  double x2 = mTofZEdge[tray-1][module-1][cell-1][iZBin+1];
1769  double y1 = mTofZCorr[tray-1][module-1][cell-1][iZBin];
1770  double y2 = mTofZCorr[tray-1][module-1][cell-1][iZBin+1];
1771  double dcorr = y1 + (z-x1)*(y2-y1)/(x2-x1);
1772 
1773  tofcorr -= dcorr;
1774  LOG_DEBUG << "zHit correction: "<<dcorr<<endm;
1775 
1776  } else {
1777  LOG_DEBUG << " Z out of range! EXIT! " << endm;
1778  return -9999.;
1779  }
1780 
1781  }
1782 
1783  LOG_DEBUG << " Corrected tof: tofcorr = " << tofcorr << endm;
1784 
1785 // Special handling of Run15 slewing correction (Bassam)
1786  const double SLEWCORR = 0.19; // slewing correction constant for Run 15 in ns
1787  if(mRun15Slew && tot<15.) {
1788  if(tot<13.) {
1789  tofcorr -= SLEWCORR;
1790  }
1791  else {
1792  tofcorr -= (SLEWCORR/4.)*(15.-tot)*(15.-tot);
1793  }
1794  }
1795 
1796  if(mRun15Slew && tot>19. && tot<26.) {
1797  double delta = (22.5-tot)/3.5;
1798  tofcorr += .110*(1.0-delta*delta);
1799  }
1800 
1801  if (mRun15Slew) LOG_DEBUG << " Corrected tof (Run15Slew): tofcorr = " << tofcorr << endm;
1802 
1803 
1804  return tofcorr;
1805 }
1806 
1807 //_____________________________________________________________________________
1808 void StBTofCalibMaker::tstart(const Double_t vz, Double_t *tstart, Double_t *tdiff)
1809 {
1810  if ( mForceTStartZero ){
1811  *tstart = 0;
1812  return;
1813  }
1814  *tstart = -9999.;
1815  *tdiff = -9999.;
1816 
1817  if(fabs(vz)>200.) {LOG_INFO << "tstart: vz too big" << endm; return;}
1818 
1819  Double_t TSum = mTSumEast + mTSumWest;
1820 
1821  if(mNEast+mNWest>0) {
1822  *tstart = (TSum-(mNEast-mNWest)*vz/(C_C_LIGHT/1.e9))/(mNEast+mNWest);
1823  }
1824  if ( mNEast>=mVPDEastHitsCut && mNWest>=mVPDWestHitsCut ) {
1825  *tdiff = (mTSumEast/mNEast - mTSumWest/mNWest)/2. - vz/(C_C_LIGHT/1.e9);
1826  }
1827 
1828  return;
1829 }
1830 
1831 //_____________________________________________________________________________
1832 void StBTofCalibMaker::tstart_NoVpd(const StBTofCollection *btofColl, const StPrimaryVertex *pVtx, Double_t *tstart)
1833 {
1834  if ( mForceTStartZero ){
1835  *tstart = 0;
1836  return;
1837  }
1838 
1839  *tstart = -9999.;
1840  if(!btofColl) return;
1841 
1842  const StSPtrVecBTofHit &tofHits = btofColl->tofHits();
1843  Int_t nCan = 0;
1844  Double_t tSum = 0.;
1845  Double_t t0[5000];
1846  Double_t tCanFirst = 99999.;
1847  Double_t tCanLast = -99999.;
1848 
1849  memset(t0, 0., sizeof(t0));
1850  for(size_t i=0;i<tofHits.size();i++) {
1851  StBTofHit *aHit = dynamic_cast<StBTofHit*>(tofHits[i]);
1852  if(!aHit) continue;
1853 
1854  isMcFlag = kFALSE; // Check to see if the hit is simulated.
1855  if ( aHit->qaTruth() == 1 ) {
1856  isMcFlag = kTRUE;
1857  }
1858 
1859  int trayId = aHit->tray();
1860  if(trayId>0&&trayId<=mNTray) {
1861  StGlobalTrack *gTrack = dynamic_cast<StGlobalTrack*>(aHit->associatedTrack());
1862  if(!gTrack) continue;
1863  StPrimaryTrack *pTrack = dynamic_cast<StPrimaryTrack*>(gTrack->node()->track(primary));
1864  if(!pTrack) continue;
1865  if(pTrack->vertex() != pVtx) continue;
1866  StThreeVectorF mom = pTrack->geometry()->momentum();
1867  double ptot = mom.mag();
1868  int q = pTrack->geometry()->charge();
1869 
1870  static StTpcDedxPidAlgorithm PidAlgorithm;
1871  static StPionPlus* Pion = StPionPlus::instance();
1872  static StProton* Proton = StProton::instance();
1873  static StKaonPlus* Kaon = StKaonPlus::instance();
1874  static StElectron* Electron = StElectron::instance();
1875 
1876  const StParticleDefinition* pd = pTrack->pidTraits(PidAlgorithm);
1877  double nSigPi = -999.;
1878  double nSigP = -999.;
1879  double nSigKaon = -999.;
1880  double nSigElectron = -999.;
1881  float nHitsDedx = 0.;
1882 
1883  if(pd) {
1884  nSigPi = PidAlgorithm.numberOfSigma(Pion);
1885  nSigP = PidAlgorithm.numberOfSigma(Proton);
1886  nSigKaon = PidAlgorithm.numberOfSigma(Kaon);
1887  nSigElectron = PidAlgorithm.numberOfSigma(Electron);
1888 
1889  nHitsDedx = PidAlgorithm.traits()->numberOfPoints();
1890  }
1891 
1892  float nHitsFit = pTrack->fitTraits().numberOfFitPoints();
1893  float nHitsPoss = pTrack->numberOfPossiblePoints();
1894 
1895  double pT = mom.perp();
1896  double eta = mom.pseudoRapidity();
1897  double zVtx = pVtx->position().z();
1898 
1899  const StPtrVecTrackPidTraits& theTofPidTraits = pTrack->pidTraits(kTofId);
1900  if(!theTofPidTraits.size()) continue;
1901 
1902  StTrackPidTraits *theSelectedTrait = theTofPidTraits[theTofPidTraits.size()-1];
1903  if(!theSelectedTrait) continue;
1904 
1905  StBTofPidTraits *pidTof = dynamic_cast<StBTofPidTraits *>(theSelectedTrait);
1906  if(!pidTof) continue;
1907 
1908  double tot = aHit->tot(); // ns
1909  double tof = aHit->leadingEdgeTime();
1910  double zhit = pidTof->zLocal();
1911 
1912  int moduleChan = (aHit->module()-1)*6 + (aHit->cell()-1);
1913  Double_t tofcorr = tof;
1914  if (!isMcFlag) {
1915  tofcorr = tofAllCorr(tof, tot, zhit, trayId, moduleChan);
1916  }
1917  if(tofcorr<0.) continue;
1918 
1919  StThreeVectorF primPos = pVtx->position();
1920  StPhysicalHelixD helix = pTrack->geometry()->helix();
1921  double L = tofPathLength(&primPos, &pidTof->position(), helix.curvature());
1922  double tofPi = L*sqrt(M_PION_PLUS*M_PION_PLUS+ptot*ptot)/(ptot*(C_C_LIGHT/1.e9));
1923  double tofP = L*sqrt(M_PROTON*M_PROTON+ptot*ptot)/(ptot*(C_C_LIGHT/1.e9));
1924 
1925  // use lose cut for low energies to improve the efficiency
1926  if(mFXTMode) //use both Pion and Proton for FXT mode
1927  {
1928  if(( (q<0 && ptot>0.2) || (q>0 && ptot>0.2 && ptot<1.0)) && fabs(nSigPi)<2.0)//pi selection
1929  {
1930  tSum += tofcorr - tofPi;
1931  t0[nCan] = tofcorr - tofPi;
1932  nCan++;
1933  }
1934  else if(q>0 && fabs(nSigP)<2.0)//proton selection
1935  {
1936  tSum += tofcorr - tofP;
1937  t0[nCan] = tofcorr - tofP;
1938  nCan++;
1939  }
1940  }
1941  else if(mPPPAMode || mPPPAPionSel)
1942  {
1943  // Obtaining gDCA
1944  const StDcaGeometry *dcaGeometry = gTrack->dcaGeometry();
1945  if(!dcaGeometry) continue;
1946  Double_t vtx[3] = {primPos[0], primPos[1], primPos[2]};
1947  THelixTrack thelix = dcaGeometry->thelix();
1948  thelix.Move(thelix.Path(vtx)); // move along the helix to the point that is closest to the vertex
1949  const Double_t *pos = thelix.Pos();
1950  StThreeVectorF gDca = StThreeVectorF(pos[0],pos[1],pos[2]) - pVtx->position();
1951  double gDcaMag = gDca.mag();
1952  if(tot > 25.) {
1953  continue;
1954  }
1955  if(ptot>0.2 && ptot<0.7 && fabs(nSigPi)< 2.0
1956  && nHitsFit>20
1957  && nHitsFit>0.51*nHitsPoss
1958  && nHitsDedx>0.5*nHitsFit
1959  && fabs(gDcaMag)<2.0
1960  && nSigKaon<-3.0
1961  && nSigElectron<-3.0
1962  && pT>0.18
1963  )
1964  {
1965  double FDGCNST = fudgeFactor(eta, zVtx);
1966  double nvrsBeta = sqrt(1.+(M_PION_PLUS/ptot)*(M_PION_PLUS/ptot));
1967  double bthBlchFrml = nvrsBeta*nvrsBeta*(BTHBLCHCNST+2.*log(ptot/M_PION_PLUS));
1968  double dTofPi = FDGCNST*L*L*((M_PION_PLUS/ptot)*(M_PION_PLUS/ptot)*(1./ptot)*bthBlchFrml);
1969  Double_t tofCorrPi = tofcorr - tofPi - dTofPi;
1970  if(tofCorrPi < tCanFirst) {
1971  tCanFirst = tofCorrPi;
1972  }
1973  if(tofCorrPi > tCanLast) {
1974  tCanLast = tofCorrPi;
1975  }
1976  tSum += tofCorrPi;
1977  t0[nCan] = tofCorrPi;
1978  nCan++;
1979  }
1980  }
1981  else//If not FXT, Only use Pion
1982  {
1983  if(ptot>0.2 && ptot<0.6 && fabs(nSigPi)< 2.0)
1984  {
1985  tSum += tofcorr - tofPi;
1986  t0[nCan] = tofcorr - tofPi;
1987  nCan++;
1988  }
1989  }
1990 
1991  }
1992  }
1993 
1994  mNTzeroCan = nCan;
1995  mTCanFirst = tCanFirst;
1996  mTCanLast = tCanLast;
1997 
1998  if(nCan<=0) {
1999  *tstart = -9999.;
2000  return;
2001  }
2002 
2003  Int_t nTzero = nCan;
2004 
2005 if(mPPPAMode || mPPPAOutlierRej || mFXTMode)
2006 {
2007  const float outlierCut = 2.5 * 0.086; // numberOfSigmas * time resolution (sigma) of a BTof pad
2008  if(nCan>1) {
2009  while(nTzero>2) {
2010  int iMax = 0;
2011  double DtiMax = 0.0;
2012  for(int i=0;i<nTzero;i++) {
2013  double tdiff = fabs(t0[i] - (tSum-t0[i])/(nTzero-1));
2014  if(tdiff>DtiMax) {
2015  iMax = i;
2016  DtiMax = tdiff;
2017  } // end of if(tdiff>DtiMax)
2018  } // end of for(int i=0;i<nTzero;i++)
2019  if(DtiMax>sqrt(1.0+1.0/(nTzero-1))*outlierCut) {
2020  tSum -= t0[iMax];
2021  t0[iMax] = t0[nTzero-1];
2022  nTzero--;
2023  } // end of if(DtiMax>1.0)
2024  else {
2025  break;
2026  }
2027  } // end of while(nTzero>2)
2028  if(nTzero==2){
2029  if(fabs(t0[0]-t0[1])>sqrt(2.0)*outlierCut) {
2030  if(mFXTMode){
2031  nTzero = 0; // in FXT mode, if only 2 tracks and too far apart in time, set nT0 to 0
2032  }else{
2033  float vpdTDiffCut = 0.25; // Acceptance window between an individual t0 vlaue
2034  // and VPD start time based on the VPD resolution
2035  double vpdStartTime = 0.;
2036  Double_t vz = pVtx->position().z();
2037  vpdTStartForBTofComparison(vz, &vpdStartTime);
2038  if(fabs(t0[0]-vpdStartTime)<vpdTDiffCut && fabs(t0[1]-vpdStartTime)<vpdTDiffCut) {
2039  if(fabs(t0[0]-vpdStartTime) < fabs(t0[1]-vpdStartTime)) {
2040  tSum -= t0[1];
2041  nTzero--;
2042  }
2043  else {
2044  tSum -= t0[0];
2045  t0[0] = t0[1];
2046  nTzero--;
2047  }
2048  } // end of if(fabs(t0[0]-vpdTStartTime)<0.25 && fabs(t0[1]-vpdTStartTime)<0.25)
2049  else if(fabs(t0[0]-vpdStartTime)<vpdTDiffCut) {
2050  tSum -= t0[1];
2051  nTzero--;
2052  }
2053  else if(fabs(fabs(t0[1]-vpdStartTime)<vpdTDiffCut)) {
2054  tSum -= t0[0];
2055  t0[0] = t0[1];
2056  nTzero--;
2057  }
2058  else {
2059  nTzero = 0;
2060  }
2061  }// end else (!mFXTMode)
2062  } // end if(fabs(t0[0]-t0[1])>sqrt(2.0)*outlierCut)
2063  } // end of if(nTzero==2)
2064  } // end of if(nCan>1)
2065  else if(mFXTMode && nCan == 1){
2066  nTzero = 0; //in FXT mode, if only 1 track set nT0 to 0 since contamination too easy
2067  }
2068 }
2069 else
2070 {
2071  if(nCan>1) { // remove hits too far from others
2072  for(int i=0;i<nCan;i++) {
2073  double tdiff = t0[i] - (tSum-t0[i])/(nTzero-1);
2074  if(fabs(tdiff)>5.0) {
2075  tSum -= t0[i];
2076  nTzero--;
2077  }
2078  }
2079  }
2080  }
2081 
2082  mNTzero = nTzero;
2083 
2084  *tstart = nTzero>0 ? tSum / nTzero : -9999.;
2085 
2086  return;
2087 }
2088 
2089 //_____________________________________________________________________________
2090 void StBTofCalibMaker::tstart_NoVpd(const StMuDst *muDst, const StMuPrimaryVertex *pVtx, Double_t *tstart)
2091 {
2092  if ( mForceTStartZero ){
2093  *tstart = 0;
2094  return;
2095  }
2096  *tstart = -9999.;
2097  if(!muDst) return;
2098 
2099  Int_t nBTofHits = muDst->numberOfBTofHit();
2100  Int_t nCan = 0;
2101  Double_t tSum = 0.;
2102  Double_t t0[5000];
2103 
2104  Double_t tCanFirst = 99999.;
2105  Double_t tCanLast = -99999.;
2106  Double_t gDCA[5000];
2107  Double_t ptotSave[5000];
2108  Double_t etaSave[5000];
2109 
2110  memset(t0, 0., sizeof(t0));
2111  int calibfailed(0); // keep counter
2112  for(int i=0;i<nBTofHits;i++) {
2113  StMuBTofHit *aHit = (StMuBTofHit*)muDst->btofHit(i);
2114  if(!aHit) {
2115  continue;
2116  }
2117 
2118  isMcFlag = kFALSE; // Check to see if the hit is simulated.
2119  if ( aHit->qaTruth() == 1 ) {
2120  isMcFlag = kTRUE;
2121  }
2122 
2123  int trayId = aHit->tray();
2124  if(trayId>0&&trayId<=mNTray) {
2125  StMuTrack *gTrack = aHit->globalTrack();
2126  if(!gTrack) continue;
2127  StMuTrack *pTrack = aHit->primaryTrack();
2128  if(!pTrack) continue;
2129  StMuPrimaryVertex *aVtx = muDst->primaryVertex(pTrack->vertexIndex());
2130  if(aVtx != pVtx) continue;
2131  StThreeVectorF mom = pTrack->momentum();
2132  double ptot = mom.mag();
2133  int q = pTrack->charge();
2134 
2135  double nSigPi = pTrack->nSigmaPion();
2136  double nSigP = pTrack->nSigmaProton();
2137 
2138  float nHitsFit = pTrack->nHitsFit();
2139  float nHitsPoss = pTrack->nHitsPoss();
2140  float nHitsDedx = pTrack->nHitsDedx();
2141 
2142  double pT = pTrack->pt();
2143  double nSigKaon = pTrack->nSigmaKaon();
2144  double nSigElectron = pTrack->nSigmaElectron();
2145  double gDcaMag = pTrack->dcaGlobal().mag();
2146  double eta = pTrack->eta();
2147  double zVtx = pVtx->position().z();
2148 
2149  StMuBTofPidTraits pidTof = pTrack->btofPidTraits();
2150 
2151  double tot = aHit->tot(); // ns
2152  double tof = aHit->leadingEdgeTime();
2153  double zhit = pidTof.zLocal();
2154 
2155  int moduleChan = (aHit->module()-1)*6 + (aHit->cell()-1);
2156 
2157  Double_t tofcorr = tof;
2158  if (!isMcFlag) {
2159  tofcorr = tofAllCorr(tof, tot, zhit, trayId, moduleChan);
2160  if(tofcorr<0.) {
2161 // LOG_WARN << " Calibration failed! ... " << endm;
2162  calibfailed++; continue;
2163  }
2164  }
2165  if(tofcorr<0.) continue;
2166 
2167  StThreeVectorF primPos = pVtx->position();
2168  StPhysicalHelixD helix = pTrack->helix();
2169  double L = tofPathLength(&primPos, &pidTof.position(), helix.curvature());
2170  double tofPi = L*sqrt(M_PION_PLUS*M_PION_PLUS+ptot*ptot)/(ptot*(C_C_LIGHT/1.e9));
2171  double tofP = L*sqrt(M_PROTON*M_PROTON+ptot*ptot)/(ptot*(C_C_LIGHT/1.e9));
2172 
2173  // For low energies, lose cut to improve the efficiency in peripheral collisions - resolution should be not a big issue
2174  if(mFXTMode) //use both Pion and Proton for FXT mode
2175  {
2176  if(( (q<0 && ptot>0.2) || (q>0 && ptot>0.2 && ptot<1.0)) && fabs(nSigPi)<2.0)//pi selection
2177  {
2178  tSum += tofcorr - tofPi;
2179  t0[nCan] = tofcorr - tofPi;
2180  nCan++;
2181  }
2182  else if(q>0 && fabs(nSigP)<2.0)//proton selection
2183  {
2184  tSum += tofcorr - tofP;
2185  t0[nCan] = tofcorr - tofP;
2186  nCan++;
2187  }
2188  }
2189  // calculating the total time, individual time,
2190  // and multiplicity of candidate particles for the pppA mode
2191  else if(mPPPAMode || mPPPAPionSel)
2192  {
2193  if(tot>25.) {
2194  continue;
2195  }
2196  if(ptot>0.2 && ptot<0.7 && fabs(nSigPi)< 2.0
2197  && nHitsFit>20
2198  && nHitsFit>0.51*nHitsPoss
2199  && nHitsDedx>0.5*nHitsFit
2200  && fabs(gDcaMag)<2.0
2201  && nSigKaon<-3.0
2202  && nSigElectron<-3.0
2203  && pT>0.18)
2204  {
2205  double FDGCNST = fudgeFactor(eta, zVtx);
2206  double nvrsBeta = sqrt(1.+(M_PION_PLUS/ptot)*(M_PION_PLUS/ptot));
2207  double bthBlchFrml = nvrsBeta*nvrsBeta*(BTHBLCHCNST+2.*log(ptot/M_PION_PLUS));
2208  double dTofPi = FDGCNST*L*L*((M_PION_PLUS/ptot)*(M_PION_PLUS/ptot)*(1./ptot)*bthBlchFrml);
2209  Double_t tofCorrPi = tofcorr - tofPi - dTofPi;
2210  if(tofCorrPi < tCanFirst) {
2211  tCanFirst = tofCorrPi;
2212  }
2213  if(tofCorrPi > tCanLast) {
2214  tCanLast = tofCorrPi;
2215  }
2216  tSum += tofCorrPi;
2217  t0[nCan] = tofCorrPi;
2218  gDCA[nCan] = gDcaMag;
2219  ptotSave[nCan] = ptot;
2220  etaSave[nCan] = eta;
2221  nCan++;
2222  }
2223  }
2224  else//If not FXT, Only use Pion
2225  {
2226  if(ptot>0.2 && ptot<0.6 && fabs(nSigPi)< 2.0)
2227  {
2228  tSum += tofcorr - tofPi;
2229  t0[nCan] = tofcorr - tofPi;
2230  nCan++;
2231  }
2232  }
2233 
2234  }//tray
2235  }//tof hit
2236 
2237  if (calibfailed) LOG_WARN << "tstart_NoVpd calibrations failures: " << calibfailed << " (out of "<< nBTofHits << "BTOF hits)"<< endm;
2238 
2239  mNTzeroCan = nCan;
2240  mTCanFirst = tCanFirst;
2241  mTCanLast = tCanLast;
2242 
2243 
2244  if(nCan<=0) {
2245  *tstart = -9999.;
2246  return;
2247  }
2248 
2249  Int_t nTzero = nCan;
2250 
2251  // Adjusting the default outlier rejection to accomodate for PPPAMode
2252  if(mPPPAMode || mPPPAOutlierRej || mFXTMode)
2253  {
2254  if(mPPPAModeHist) {
2255  if(nCan==1) { // filling global DCA histograms before outlier rejection
2256  hGDcaArray[0]->Fill(fabs(gDCA[0]));
2257  }
2258  else if(nCan==2) {
2259  for(int i=0;i<nCan;i++) {
2260  hGDcaArray[1]->Fill(fabs(gDCA[i]));
2261  }
2262  }
2263  else if(nCan>2) {
2264  for(int i=0;i<nCan;i++) {
2265  hGDcaArray[2]->Fill(fabs(gDCA[i]));
2266  }
2267  }
2268  }
2269  const float outlierCut = 2.5 * 0.086; // numberOfSigmas * time resolution (sigma) of a BTof pad
2270  if(nCan>1) {
2271  while(nTzero>2) {
2272  int iMax = 0;
2273  double DtiMax = 0.0;
2274  for(int i=0;i<nTzero;i++) {
2275  double tdiff = fabs(t0[i] - (tSum-t0[i])/(nTzero-1));
2276  if(tdiff>DtiMax) {
2277  iMax = i;
2278  DtiMax = tdiff;
2279  } // end of if(tdiff>DtiMax)
2280  } // end of for(int i=0;i<nTzero;i++)
2281  if(DtiMax>sqrt(1.0+1.0/(nTzero-1))*outlierCut) {
2282  tSum -= t0[iMax];
2283  t0[iMax] = t0[nTzero-1];
2284  gDCA[iMax] = gDCA[nTzero-1];
2285  ptotSave[iMax] = ptotSave[nTzero-1];
2286  etaSave[iMax] = etaSave[nTzero-1];
2287  nTzero--;
2288  } // end of if(DtiMax>1.0)
2289  else {
2290  break;
2291  }
2292  } // end of while(nTzero>2)
2293  if(nTzero==2){
2294  if(fabs(t0[0]-t0[1])>sqrt(2.0)*outlierCut) {
2295  if(mFXTMode){
2296  nTzero = 0;
2297  }else{
2298  float vpdTDiffCut = 0.25; // Acceptance window between an individual t0 vlaue
2299  //and VPD start time based on the VPD resolution
2300  double vpdStartTime = 0.;
2301  Double_t vz = pVtx->position().z();
2302  vpdTStartForBTofComparison(vz, &vpdStartTime);
2303  if(fabs(t0[0]-vpdStartTime)<vpdTDiffCut && fabs(t0[1]-vpdStartTime)<vpdTDiffCut) {
2304  if(fabs(t0[0]-vpdStartTime) < fabs(t0[1]-vpdStartTime)) {
2305  tSum -= t0[1];
2306  nTzero--;
2307  }
2308  else {
2309  tSum -= t0[0];
2310  t0[0] = t0[1];
2311  gDCA[0] = gDCA[1];
2312  ptotSave[0] = ptotSave[1];
2313  etaSave[0] = etaSave[1];
2314  nTzero--;
2315  }
2316  } // end of if(fabs(t0[0]-vpdTStartTime)<0.25 && fabs(t0[1]-vpdTStartTime)<0.25)
2317  else if(fabs(t0[0]-vpdStartTime)<vpdTDiffCut) {
2318  tSum -= t0[1];
2319  nTzero--;
2320  }
2321  else if(fabs(fabs(t0[1]-vpdStartTime)<vpdTDiffCut)) {
2322  tSum -= t0[0];
2323  t0[0] = t0[1];
2324  gDCA[0] = gDCA[1];
2325  ptotSave[0] = ptotSave[1];
2326  etaSave[0] = etaSave[1];
2327  nTzero--;
2328  }
2329  else {
2330  nTzero = 0;
2331  }
2332  }// end else (!mFXTMode)
2333  }//end if(fabs(t0[0]-t0[1])>sqrt(2.0)*outlierCut)
2334  } // end of if(nTzero==2)
2335  } // end of if(nCan>1)
2336  else if(mFXTMode && nCan == 1){
2337  nTzero = 0; //in FXT mode, if only 1 track set nT0 to 0 since contamination too easy
2338  }
2339  }
2340  else
2341  {
2342  if(nCan>1) { // remove hits too far from others
2343  for(int i=0;i<nCan;i++) {
2344  double tdiff = t0[i] - (tSum-t0[i])/(nTzero-1);
2345  if(fabs(tdiff)>5.0) {
2346  tSum -= t0[i];
2347  nTzero--;
2348  }
2349  }
2350  }
2351  }
2352 
2353  mNTzero = nTzero;
2354 
2355  *tstart = nTzero>0 ? tSum / nTzero : -9999.;
2356 
2357  if(mPPPAModeHist) {
2358  if(nTzero>1) {
2359  for(int i=0;i<nTzero; i++) {
2360  double Dti = t0[i] - (tSum-t0[i])/(nTzero-1);
2361  if(nTzero>10) {
2362  hDtiArray[9]->Fill(Dti);
2363  }
2364  else {
2365  hDtiArray[nTzero-2]->Fill(Dti);
2366  }
2367  //if(nTzero>4 && fabs(pVtx->position().perp())<1. && fabs(pVtx->position().z())<75.) {
2368  // hDtiVsPtot[0]->Fill(ptotSave[i], Dti);
2369  //}
2370  //if(nTzero>4 && fabs(pVtx->position().perp())<1. && fabs(pVtx->position().z())>75.) {
2371  // hDtiVsPtot[1]->Fill(ptotSave[i], Dti);
2372  //}
2373  double rVtx = pVtx->position().perp();
2374  double zVtx = pVtx->position().z();
2375  if(nTzero>4 && (fabs(rVtx)<1. || (fabs(zVtx-200.)<3. && fabs(rVtx)<3.))) {
2376  fillDtiVsPtotProfiles(etaSave[i], zVtx, ptotSave[i], Dti);
2377  }
2378  } // end of for(int i=0;i<mNTzero; i++)
2379  } // end of if(mNTzero>1)
2380  if(nCan==2) { // filling global DCA histograms after outlier rejection
2381  for(int i=0;i<nTzero;i++) {
2382  hGDcaArray[3]->Fill(fabs(gDCA[i]));
2383  }
2384  }
2385  else if(nCan>2) {
2386  for(int i=0;i<nTzero;i++) {
2387  hGDcaArray[4]->Fill(fabs(gDCA[i]));
2388  }
2389  }
2390  } // end of if(mPPPAModeHist)
2391 
2392 
2393  return;
2394 }
2395 
2396 //_____________________________________________________________________________
2397 void StBTofCalibMaker::vpdTStartForBTofComparison(Double_t vz, Double_t *vpdStartTime)
2398 {
2399  if(!mBTofHeader){
2400  *vpdStartTime = -9999.0;
2401  return;
2402  }
2403  double vzNs = vz/mC_Light;
2404  //int nVpdWCandHits = 0;
2405  //int nVpdECandHits = 0;
2406  Double_t vpdHitTStart[38]; // Array will contain the start time of West + East VPD hit start time
2407  //Bool_t isEastHit[38]; // true if the hit came from VPDE false if hit came from VPDW
2408 
2409  int nVpdTotCandHits = 0; // Total number of West + East VPD candidate hits in a given event
2410  double tSum = 0.;
2411 
2412  for(int i=0; i<19; i++) {
2413  double vpdTime = mBTofHeader->vpdTime(west,i+1); // Here, the index has to be between 1-19
2414  //since the btofHeader function defines tubeId as tubeId-1.
2415  if(vpdTime<1.e-4) {
2416  continue;
2417  }
2418  vpdHitTStart[nVpdTotCandHits] = vpdTime + vzNs - DEDXTCORR[iYr];
2419  //isEastHit[nVpdTotCandHits] = kFALSE;
2420  tSum += vpdHitTStart[nVpdTotCandHits];
2421  //nVpdWCandHits++;
2422  nVpdTotCandHits++;
2423  }
2424 
2425  for(int i=0; i<19; i++) {
2426  double vpdTime = mBTofHeader->vpdTime(east,i+1); // Here, the index has to be between 1-19
2427  //since the btofHeader function defines tubeId as tubeId-1.
2428  if(vpdTime<1.e-4) {
2429  continue;
2430  }
2431  vpdHitTStart[nVpdTotCandHits] = vpdTime - vzNs - DEDXTCORR[iYr];
2432  //isEastHit[nVpdTotCandHits] = kTRUE;
2433  tSum += vpdHitTStart[nVpdTotCandHits];
2434  //nVpdECandHits++;
2435  nVpdTotCandHits++;
2436  }
2437 
2438  if(nVpdTotCandHits<1) {
2439  //*mNVpdWHits = 0;
2440  //*mNVpdEHits = 0;
2441  //*mNVpdTzero = 0;
2442  *vpdStartTime = -9999.0;
2443  return;
2444  }
2445 
2446  Int_t nVpdTzero = nVpdTotCandHits;
2447 
2448  const float vpdTStartOulierCut = 2.5*0.130; // numberOfSigmas * time resolution of a vpd phototube
2449 
2450  if(nVpdTotCandHits>1) {
2451  while(nVpdTzero>2) {
2452  int iVpdMax = 0;
2453  double vpdDtiMax = 0.;
2454  for(int i=0; i<nVpdTzero; i++) {
2455  double vpdTdiff = fabs(vpdHitTStart[i] - (tSum-vpdHitTStart[i])/(nVpdTzero-1));
2456  if(vpdTdiff>vpdDtiMax) {
2457  iVpdMax = i;
2458  vpdDtiMax = vpdTdiff;
2459  } // end of if(vpdTdiff>vpdDtiMax)
2460  } // end of for(int i=0; i<nVpdTzero; i++)
2461  if(vpdDtiMax > sqrt(1.0+1.0/(nVpdTzero-1))*vpdTStartOulierCut) {
2462  tSum -= vpdHitTStart[iVpdMax];
2463  vpdHitTStart[iVpdMax] = vpdHitTStart[nVpdTzero-1];
2464  //isEastHit[iVpdMax] = isEastHit[nVpdTzero-1];
2465  nVpdTzero--;
2466  } // end of if(vpdDitMax>sqrt(1.0+1.0/(nVpdTzero-1))*vpdTStartOulierCut)
2467  else{
2468  break;
2469  }
2470  } // end of while(nVpdTzero)
2471  if(nVpdTzero==2) {
2472  if(fabs(vpdHitTStart[0]-vpdHitTStart[1]) > sqrt(2.0)*vpdTStartOulierCut) {
2473  nVpdTzero = 0;
2474  }
2475  } // end of if(nTzero==2)
2476  } // end of nVpdTzero
2477 
2478  *vpdStartTime = nVpdTzero>0 ? tSum/nVpdTzero : -9999.;
2479 }
2480 
2481 //_____________________________________________________________________________
2482 double StBTofCalibMaker::fudgeFactor(double eta, double zVtx)
2483 {
2484  double ffArr[2][22] = {
2485  {6.1e-8, 6.95e-8, 7.9e-8, 11.2e-8, 14.7e-8, 13.2e-8, 10.7e-8, 8.8e-8, 8.21e-8, 8.31e-8, 9.01e-8, 8.31e-8, 7.2e-8, 9.e-8, 8.8e-8, 12.e-8, 8.9e-8, 6.e-8, 5.1e-8, 4.7e-8, 5.9e-8, 3.3e-8},
2486  {4.0e-8, 3.25e-8, 4.45e-8, 4.75e-8, 4.1e-8, 3.6e-8, 3.8e-8, 3.7e-8, 3.6e-8, 3.5e-8, 3.8e-8, 3.7e-8, 3.5e-8, 3.8e-8, 3.3e-8, 4.0e-8, 4.4e-8, 4.5e-8, 2.7e-8, 3.5e-8, 4.83e-8, 3.3e-8}
2487  };
2488  double FF = 0.;
2489  if(eta<=0.) {
2490  if(zVtx<-75.) {
2491  FF = ffArr[iYr][0];
2492  }
2493  else if(zVtx<-50.) {
2494  FF = ffArr[iYr][2];
2495  }
2496  else if(zVtx<-30.) {
2497  FF = ffArr[iYr][4];
2498  }
2499  else if(zVtx<-10.) {
2500  FF = ffArr[iYr][6];
2501  }
2502  else if(zVtx<0.) {
2503  FF = ffArr[iYr][8];
2504  }
2505  else if(zVtx<10.) {
2506  FF = ffArr[iYr][10];
2507  }
2508  else if(zVtx<30.) {
2509  FF = ffArr[iYr][12];
2510  }
2511  else if(zVtx<50.) {
2512  FF = ffArr[iYr][14];
2513  }
2514  else if(zVtx<75.) {
2515  FF = ffArr[iYr][16];
2516  }
2517  else if(zVtx<196.){
2518  FF = ffArr[iYr][18];
2519  }
2520  else {
2521  if (eta < -1.) {
2522  FF = ffArr[iYr][20];
2523  }
2524  else {
2525  FF = ffArr[iYr][21];
2526  }
2527  }
2528  } // end of if(eta<=0.)
2529  else {
2530  if(zVtx<-75.) {
2531  FF = ffArr[iYr][1];
2532  }
2533  else if(zVtx<-50.) {
2534  FF = ffArr[iYr][3];
2535  }
2536  else if(zVtx<-30.) {
2537  FF = ffArr[iYr][5];
2538  }
2539  else if(zVtx<-10.) {
2540  FF = ffArr[iYr][7];
2541  }
2542  else if(zVtx<0.) {
2543  FF = ffArr[iYr][9];
2544  }
2545  else if(zVtx<10.) {
2546  FF = ffArr[iYr][11];
2547  }
2548  else if(zVtx<30.) {
2549  FF = ffArr[iYr][13];
2550  }
2551  else if(zVtx<50.) {
2552  FF = ffArr[iYr][15];
2553  }
2554  else if(zVtx<75.) {
2555  FF = ffArr[iYr][17];
2556  }
2557  else if(zVtx<196.){
2558  FF = ffArr[iYr][19];
2559  }
2560  else {
2561  FF = ffArr[iYr][21];
2562  }
2563  }
2564  return FF;
2565 }
2566 
2567 //_____________________________________________________________________________
2568 void StBTofCalibMaker::fillDtiVsPtotProfiles(double eta, double zVtx, double ptot, double Dti)
2569 {
2570  int etaIdx;
2571  int zIdx;
2572 
2573  if(zVtx<-75.) {
2574  zIdx = 0;
2575  }
2576  else if(zVtx<-50.) {
2577  zIdx = 1;
2578  }
2579  else if(zVtx<-30.) {
2580  zIdx = 2;
2581  }
2582  else if(zVtx<-10.) {
2583  zIdx = 3;
2584  }
2585  else if(zVtx<0.) {
2586  zIdx = 4;
2587  }
2588  else if(zVtx<10.) {
2589  zIdx = 5;
2590  }
2591  else if(zVtx<30.) {
2592  zIdx = 6;
2593  }
2594  else if(zVtx<50.) {
2595  zIdx = 7;
2596  }
2597  else if(zVtx<75.) {
2598  zIdx = 8;
2599  }
2600  else if(fabs(zVtx-200.)<3.){
2601  zIdx = 10;
2602  }
2603  else {
2604  zIdx = 9;
2605  }
2606 
2607  if(eta<=0.) {
2608  if(zIdx==10) {
2609  if(eta<-1.) {
2610  etaIdx = 0;
2611  }
2612  else {
2613  etaIdx = 1;
2614  }
2615  }
2616  else {
2617  etaIdx = 0;
2618  }
2619  }
2620  else {
2621  etaIdx = 1;
2622  }
2623 
2624  hDtiVsPtot[zIdx][etaIdx]->Fill(ptot, Dti);
2625 }
2626 
2627 //_____________________________________________________________________________
2628 void StBTofCalibMaker::archiveVpdHitInfo()
2629 {
2630  if(!mBTofHeader) return;
2631  if(mEvtVtxZ<-999.) return;
2632 
2633  const double SIGMA = 0.13; // ns
2634  float dist = 2.5;
2635 
2636  double vzNs = mEvtVtxZ/mC_Light; // mC_Light is c in cm/ns
2637 
2638  int wHitCan = 0;
2639  int eHitCan = 0;
2640  int goodWHits = 0;
2641  int goodEHits = 0;
2642 
2643  double wEarlHit = 99999.;
2644  double eEarlHit = 99999.;
2645 
2646  double wClosestHit = 99999.;
2647  double wClosestDiff = 99999.;
2648  double eClosestHit = 99999.;
2649  double eClosestDiff = 99999.;
2650 
2651  double wLateHit = -99999.;
2652  double eLateHit = -99999.;
2653 
2654  for(int i=0; i<mNVPD; i++) {
2655  double wHitTime = mBTofHeader->vpdTime(west,i+1);
2656  if(wHitTime<1.e-4) {
2657  continue;
2658  }
2659  wHitCan++;
2660  double wHitTimeCorr = wHitTime + vzNs - DEDXTCORR[iYr];
2661  if(wHitTimeCorr<wEarlHit) {
2662  wEarlHit = wHitTimeCorr;
2663  }
2664  if(wHitTimeCorr>wLateHit) {
2665  wLateHit = wHitTimeCorr;
2666  }
2667  if(mTStart>-1000.) {
2668  double wHitTDiff = wHitTimeCorr - mTStart;
2669  if(fabs(wHitTDiff)<dist*SIGMA) {
2670  goodWHits++;
2671  }
2672  if(fabs(wHitTDiff)<fabs(wClosestDiff)) {
2673  wClosestDiff = wHitTDiff;
2674  wClosestHit = wHitTimeCorr;
2675  }
2676  } // End of if(mTStart>-1000.)
2677  } // end of for(int i=0; i<MAXVPD; i++)
2678 
2679  for(int i=0; i<mNVPD; i++) {
2680  double eHitTime = mBTofHeader->vpdTime(east,i+1);
2681  if(eHitTime<1.e-4) {
2682  continue;
2683  }
2684  eHitCan++;
2685  double eHitTimeCorr = eHitTime - vzNs - DEDXTCORR[iYr];
2686  if(eHitTimeCorr<eEarlHit) {
2687  eEarlHit = eHitTimeCorr;
2688  }
2689  if(eHitTimeCorr>eLateHit) {
2690  eLateHit = eHitTimeCorr;
2691  }
2692  if(mTStart>-1000.) {
2693  double eHitTDiff = eHitTimeCorr - mTStart;
2694  if(fabs(eHitTDiff)<dist*SIGMA) {
2695  goodEHits++;
2696  }
2697  if(fabs(eHitTDiff)<fabs(eClosestDiff)) {
2698  eClosestDiff = eHitTDiff;
2699  eClosestHit = eHitTimeCorr;
2700  }
2701  } // End of if(mTStart>-1000.)
2702  } // End of for(int i=0; i<MAXVPD; i++)
2703 
2704  mEarliestVpdEHit = eEarlHit;
2705  mClosestVpdEHit = eClosestHit;
2706  mLatestVpdEHit = eLateHit;
2707 
2708  mEarliestVpdWHit = wEarlHit;
2709  mClosestVpdWHit = wClosestHit;
2710  mLatestVpdWHit = wLateHit;
2711 
2712  mVpdEGoodHits = goodEHits;
2713  mVpdWGoodHits = goodWHits;
2714  mVpdEHits = eHitCan;
2715  mVpdWHits = wHitCan;
2716 
2717  return;
2718 }
2719 
2720 //_____________________________________________________________________________
2721 void StBTofCalibMaker::bookHistograms()
2722 {
2723  hEventCounter = new TH1D("eventCounter","eventCounter",20,0,20);
2724 }
2725 
2726 //_____________________________________________________________________________
2727 void StBTofCalibMaker::writeHistograms()
2728 {
2729  // Output file
2730  TFile *theHistoFile = new TFile(mHistoFileName.c_str(), "RECREATE");
2731  LOG_INFO << "StBTofCalibMaker::writeHistograms()"
2732  << " histogram file " << mHistoFileName << endm;
2733 
2734  theHistoFile->cd();
2735 
2736  if(mHisto) {
2737  hEventCounter->Write();
2738  }
2739  return;
2740 }
2741 
2742 //_____________________________________________________________________________
2743 void StBTofCalibMaker::bookPPPAHistograms()
2744 {
2745  for(int i=0;i<10;i++) {
2746  ostringstream histogramName;
2747  ostringstream histogramTitle;
2748 
2749  if(i == 9)
2750  {
2751  histogramName << "DtiDistFornTzeroGreater10";
2752  histogramTitle << "Dti Dist. For nTzero > 10";
2753  }
2754  else
2755  {
2756  histogramName << "DtiDistFornTzero" << i+2;
2757  histogramTitle << "Dti Dist. for nTzero " << i+2;
2758  }
2759  hDtiArray[i] = new TH1D(histogramName.str().c_str(), histogramTitle.str().c_str(), 200, -1.0, 1.0);
2760  } // end of for(int i=0;i<5;i++)
2761 
2762  for(int i=0;i<5;i++) {
2763  ostringstream histogramName;
2764  ostringstream histogramTitle;
2765 
2766  if(i<2) {
2767  histogramName << "gDCABeforeOutlierRejectionForNCan" << i+1;
2768  histogramTitle << "|gDCA| Before Outlier Rejection For nCan = " << i+1;
2769  }
2770  else if(i==2) {
2771  histogramName << "gDCABeforeOutlierRejectionForNCanGreater2";
2772  histogramTitle << "|gDCA| Before Outlier Rejection For nCan > 2";
2773  }
2774  else if(i==3) {
2775  histogramName << "gDCAAfterOutlierRejectionForNCan2";
2776  histogramTitle << "|gDCA| After Outlier Rejection For nCan = 2";
2777  }
2778  else if(i==4) {
2779  histogramName << "gDCAAfterOutlierRejectionForNCanGreater2";
2780  histogramTitle << "|gDCA| After Outlier Rejection Ror nCan > 2";
2781  }
2782  hGDcaArray[i] = new TH1D(histogramName.str().c_str(), histogramTitle.str().c_str(), 200, 0., 2.5);
2783  }
2784 
2785  string zVertLocation[11] = {"z_{vtx}#LT-75", "-75#leqz_{vtx}#LT-50", "-50#leqz_{vtx}#LT-30", "-30#leqz_{vtx}#LT-10", "-10#leqz_{vtx}#LT0", "0#leqz_{vtx}#LT10", "10#leqz_{vtx}#LT30", "30#leqz_{vtx}#LT50", "50#leqz_{vtx}#LT75", "z_{vtx}#geq75 (excluding |z_{vtx}-200|<3)", "|z_{vtx}-200|<3"};
2786  string trkEta[4] = {"#eta#leq0","#eta>0", "#eta#LT-1", "#eta#geq-1"};
2787  for(int i=0; i<11; i++) {
2788  for(int j=0; j<2; j++) {
2789  ostringstream histogramName;
2790  ostringstream histogramTitle;
2791 
2792  if(i==10 && j==0) {
2793  histogramName << "hDtiVsPtot_" << i << j;
2794  histogramTitle << "Dti vs. pTot " << zVertLocation[i] << " " << trkEta[2];
2795  }
2796  else if(i==10 && j==1) {
2797  histogramName << "hDtiVsPtot_" << i << j;
2798  histogramTitle << "Dti vs. pTot " << zVertLocation[i] << " " << trkEta[3];
2799  }
2800  else {
2801  histogramName << "hDtiVsPtot_" << i << j;
2802  histogramTitle << "Dti vs. pTot " << zVertLocation[i] << " " << trkEta[j];
2803  }
2804 
2805  hDtiVsPtot[i][j] = new TProfile(histogramName.str().c_str(), histogramTitle.str().c_str(), 50, 0.2, 0.7, -1.0, 1.0);
2806  }
2807  }
2808 }
2809 
2810 //_____________________________________________________________________________
2811 void StBTofCalibMaker::writePPPAHistograms()
2812 {
2813  // Output file
2814  TFile *thePPPAHistoFile = new TFile(mPPPAModeHistoFileName.c_str(), "RECREATE");
2815 
2816  thePPPAHistoFile->cd();
2817 
2818  if(mPPPAModeHist) {
2819  for(int i=0;i<10;i++) {
2820  hDtiArray[i]->Write();
2821  }
2822  for(int i=0;i<5;i++) {
2823  hGDcaArray[i]->Write();
2824  }
2825  for(int i=0; i<11; i++) {
2826  for(int j=0; j<2; j++) {
2827  hDtiVsPtot[i][j]->Write();
2828  }
2829  }
2830  }
2831 
2832  return;
2833 }
2834 
2835 //_____________________________________________________________________________
2836 float StBTofCalibMaker::tofCellResolution(const Int_t itray, const Int_t iModuleChan)
2837 {
2838  float resolution(0.013); // 0.013 by default - 1/beta resolution
2839  if (itray<0){return resolution;}
2840 
2841  int module = iModuleChan/6 + 1;
2842  int cell = iModuleChan%6 + 1;
2843  // mBTofRes::timeres_tof() reports in picoseconds
2844  float stop_resolution = mBTofRes->timeres_tof(itray, module, cell)/1000.;
2845 
2846  float start_resolution = 0.0;
2847  if (mUseVpdStart)
2848  start_resolution = mVpdResConfig->singleTubeRes(mVPDHitPatternEast, mVPDHitPatternWest)/1000.;
2849  else
2850  start_resolution = mBTofRes->average_timeres_tof()/sqrt(mNTzero)/1000.;
2851  resolution = sqrt(stop_resolution*stop_resolution + start_resolution*start_resolution);
2852 
2853  return resolution;
2854 }
StMaker::m_Mode
Int_t m_Mode
counters
Definition: StMaker.h:81
StMuDst::primaryVertex
static StMuPrimaryVertex * primaryVertex()
return pointer to current primary vertex
Definition: StMuDst.h:322
StBTofCalibMaker::setVPDHitsCut
void setVPDHitsCut(const Int_t, const Int_t)
switch to set the VPD # of hits cut
Definition: StBTofCalibMaker.h:357
StMuDst::btofHit
static StMuBTofHit * btofHit(int i)
returns pointer to the i-th muBTofHit
Definition: StMuDst.h:419
StMuDst::globalTracks
static TObjArray * globalTracks()
returns pointer to the global tracks list
Definition: StMuDst.h:303
StBTofSimResParams::timeres_tof
double timeres_tof(unsigned int itray, unsigned int imodule, unsigned int icell)
Definition: StBTofSimResParams.h:35
StMuTrack::vertexIndex
Int_t vertexIndex() const
Returns index of associated primary vertex.
Definition: StMuTrack.cxx:600
StBTofCalibMaker::InitRun
virtual Int_t InitRun(int)
Definition: StBTofCalibMaker.cxx:322
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::pt
Double_t pt() const
Returns pT at point of dca to primary vertex.
Definition: StMuTrack.h:256
StMuTrack::nHitsDedx
UShort_t nHitsDedx() const
Return number of hits used for dEdx.
Definition: StMuTrack.h:238
StBTofCalibMaker::setHistoFileName
void setHistoFileName(const Char_t *)
set histogram output file name
Definition: StBTofCalibMaker.h:362
StVpdCalibMaker::useVpdStart
Bool_t useVpdStart() const
function for tofCalibMaker to know whether to use VPD as the start or not
Definition: StVpdCalibMaker.h:200
tof
Definition: tof.h:15
StBTofCalibMaker::setOuterGeometry
void setOuterGeometry(const bool val=kTRUE)
switch to select the Helix Geometry
Definition: StBTofCalibMaker.h:358
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
StBTofSimResParams::loadParams
void loadParams(const int runNumber=20076002)
Loads BTOF Sim Params from database.
Definition: StBTofSimResParams.h:46
StMuTrack::charge
Short_t charge() const
Returns charge.
Definition: StMuTrack.h:255
StBTofCalibMaker::Finish
virtual Int_t Finish()
Definition: StBTofCalibMaker.cxx:970
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::eta
Double_t eta() const
Returns pseudo rapidity at point of dca to primary vertex.
Definition: StMuTrack.h:257
StMuDst::primaryTracks
static TObjArray * primaryTracks()
returns pointer to a list of tracks belonging to the selected primary vertex
Definition: StMuDst.h:301
StBTofCalibMaker::StBTofCalibMaker
StBTofCalibMaker(const char *name="btofCalib")
Default constructor.
Definition: StBTofCalibMaker.cxx:178
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
StMuTrack::nHitsPoss
UShort_t nHitsPoss() const
Return number of possible hits on track.
Definition: StMuTrack.cxx:242
StMuTrack::momentum
const StThreeVectorF & momentum() const
Returns 3-momentum at dca to primary vertex.
Definition: StMuTrack.h:260
StBTofCalibMaker::~StBTofCalibMaker
virtual ~StBTofCalibMaker()
Destructor.
Definition: StBTofCalibMaker.cxx:226
kStOK
Definition: Stypes.h:40
StBTofCalibMaker::Make
virtual Int_t Make()
Definition: StBTofCalibMaker.cxx:978
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
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
THelixTrack::Move
double Move(double step)
Move along helix.
Definition: THelixTrack.cxx:452
kStErr
Definition: Stypes.h:44
StVpdSimConfig::loadVpdSimParams
void loadVpdSimParams(const int date=20160913, const int time=175725, const char *Default_time="2016-09-13 17:57:25")
Loads Vpd Sim Params from database.
Definition: StVpdSimConfig.h:109
StVpdSimConfig::singleTubeRes
double singleTubeRes(UInt_t mVPDHitPatternEast, UInt_t mVPDHitPatternWest)
Calculate correct resolution based on those tubes that were used.
Definition: StVpdSimConfig.h:44
StBTofCalibMaker::setCreateHistoFlag
void setCreateHistoFlag(Bool_t histos=kTRUE)
enable QA histogram filling
Definition: StBTofCalibMaker.h:363
TDataSet::Find
virtual TDataSet * Find(const char *path) const
Definition: TDataSet.cxx:362
Generated by   doxygen 1.8.5