StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StMtdCalibMaker.cxx
1 /*******************************************************************
2  *
3  * StMtdCalibMaker.cxx,v 1.0 2014/05/27
4  *
5  * Author: Xinjie Huang
6  *****************************************************************
7  *
8  * Description: - Mtd Calibration Maker to do the calibration for Mtd
9  *
10  * - store into StMtdPidTraits
11  *
12  *************************************************************************/
13 
14 #include <iostream>
15 #include "TH1F.h"
16 #include "TH2F.h"
17 
18 #include "StEvent.h"
19 #include "StBTofCollection.h"
20 #include "StBTofHeader.h"
21 #include "StEventTypes.h"
22 #include "Stypes.h"
23 #include "StThreeVectorD.hh"
24 #include "StHelix.hh"
25 #include "StTrackGeometry.h"
26 #include "StTrackPidTraits.h"
27 #include "PhysicalConstants.h"
28 
29 #include "StBTofCollection.h"
30 #include "StBTofHeader.h"
31 #include "StMtdCollection.h"
32 #include "StMtdHeader.h"
33 
34 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
35 #include "StMuDSTMaker/COMMON/StMuDst.h"
36 #include "StMuDSTMaker/COMMON/StMuTrack.h"
37 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
38 #include "StMuDSTMaker/COMMON/StMuEvent.h"
39 #include "StMuDSTMaker/COMMON/StMuMtdCollection.h"
40 #include "StMuDSTMaker/COMMON/StMuMtdRawHit.h"
41 #include "StMuDSTMaker/COMMON/StMuMtdHeader.h"
42 #include "StMuDSTMaker/COMMON/StMuMtdHit.h"
43 
44 #include "tables/St_mtdSlewingCorr_Table.h"
45 #include "tables/St_mtdT0Offset_Table.h"
46 #include "tables/St_mtdTriggerTimeCut_Table.h"
47 #include "tables/St_mtdPositionCorr_Table.h"
48 
49 #include "StMtdCalibMaker.h"
50 
51 ClassImp(StMtdCalibMaker)
52 
53 //_____________________________________________________________________________
54 StMtdCalibMaker::StMtdCalibMaker(const char *name) : StMaker(name)
55 {
57 
58  mStEvent = 0x0;
59  mMuDst = 0x0;
60  mDebug = kFALSE;
61  mHisto = kFALSE;
62  mUseTriggercut = kFALSE;
63  mInitFromFile = kFALSE; // default initialization from database
64  mCalibFileTot = "";
65  mCalibFileT0 = "";
66  mCalibFileTrigger = "";
67  mCalibFileDy = "";
68  mCalibFileDz = "";
69  hTimeOfFlightCorr = 0x0;
70  hAllCorr = 0x0;
71  hTimeOfFlightModule = 0x0;
72  hTimeOfFlightCorrModule = 0x0;
73  hTriggerTimeBL = 0x0;
74  hVertexzVsTpcz = 0x0;
75  hTOFTimeOfFlightTray = 0x0;
76  hDyModule = 0x0;
77  hDzModule = 0x0;
78 }
79 
80 //_____________________________________________________________________________
81 StMtdCalibMaker::~StMtdCalibMaker()
82 {
84  if(hTimeOfFlightCorr) delete hTimeOfFlightCorr;
85  if(hAllCorr) delete hAllCorr;
86  if(hTimeOfFlightModule) delete hTimeOfFlightModule;
87  if(hTimeOfFlightCorrModule) delete hTimeOfFlightCorrModule;
88  if(hTriggerTimeBL) delete hTriggerTimeBL;
89  if(hVertexzVsTpcz) delete hVertexzVsTpcz;
90  if(hTOFTimeOfFlightTray) delete hTOFTimeOfFlightTray;
91  if(hDyModule) delete hDyModule;
92  if(hDzModule) delete hDzModule;
93 }
94 
95 //____________________________________________________________________________
96 Int_t StMtdCalibMaker::Init()
97 {
98  if (mHisto) bookHistograms();
99 
100  return kStOK;
101 }
102 
103 //____________________________________________________________________________
104 Int_t StMtdCalibMaker::InitRun(Int_t runnumber)
105 {
106  //return kStOK;
107  //initialize correction table
108  Int_t dblSize = sizeof(Double_t);
109  memset(mMtdT0Corr, 0, mNBackleg*mNModule*mNCell*dblSize);
110  memset(mMtdTotCorr, 0, mNBackleg*mNModule*mNBinMax*dblSize);
111  memset(mMtdTotEdge, 0, mNBackleg*mNModule*mNBinMax*dblSize);
112  memset(mMtdDyCorr, 0, mNBackleg*mNModule*dblSize);
113  memset(mMtdDzCorr, 0, mNBackleg*mNModule*mNCell*dblSize);
114  memset(mTriggerHighEdge,0, mNBackleg*mNModule*dblSize);
115  memset(mTriggerLowEdge, 0, mNBackleg*mNModule*dblSize);
116 
117  if(!mInitFromFile) //init from database
118  {
119  LOG_INFO << "Initializing calibration parameters from data base" << endm;
120  TDataSet *mDbDataSet = 0x0;
121 
122  // extract T0 offset parameters
123  mDbDataSet = GetDataBase("Calibrations/mtd/mtdT0Offset");
124  St_mtdT0Offset* mtdT0Offset = static_cast<St_mtdT0Offset*>(mDbDataSet->Find("mtdT0Offset"));
125  if(!mtdT0Offset)
126  {
127  LOG_ERROR << "Unable to get the T0 offset parameters from data base" << endm;
128  return kStErr;
129  }
130  mtdT0Offset_st* tZero = static_cast<mtdT0Offset_st*>(mtdT0Offset->GetTable());
131  for(int j=0; j<mNBackleg*mNModule*mNCell; j++)
132  {
133  Int_t backleg = j/60;
134  Int_t module = (j-backleg*60)/12;
135  Int_t cell = j-backleg*60-module*12;
136  mMtdT0Corr[backleg][module][cell] = tZero->t0Offset[j];
137  }
138 
139  // extract Dy and Dz correction parameters
140  mDbDataSet = GetDataBase("Calibrations/mtd/mtdPositionCorr");
141  St_mtdPositionCorr* mtdPositionCorr = static_cast<St_mtdPositionCorr*>(mDbDataSet->Find("mtdPositionCorr"));
142  if(!mtdPositionCorr)
143  {
144  LOG_ERROR << "Unable to get the posistion correction parameters from data base" << endm;
145  return kStErr;
146  }
147  mtdPositionCorr_st* posCorr = static_cast<mtdPositionCorr_st*>(mtdPositionCorr->GetTable());
148  for(int j=0; j<mNBackleg*mNModule; j++)
149  {
150  Int_t backleg = j/5;
151  Int_t module = j%5;
152  mMtdDyCorr[backleg][module] = posCorr->yCorr[j];
153  }
154 
155  for(int j=0; j<mNBackleg*mNModule*mNCell; j++)
156  {
157  Int_t backleg = j/60;
158  Int_t module = (j-backleg*60)/12;
159  Int_t cell = j-backleg*60-module*12;
160  mMtdDzCorr[backleg][module][cell] = posCorr->zCorr[j];
161  }
162 
163  // extract slewing correction parameters
164  mDbDataSet = GetDataBase("Calibrations/mtd/mtdSlewingCorr");
165  St_mtdSlewingCorr* mtdSlewingCorr = static_cast<St_mtdSlewingCorr*>(mDbDataSet->Find("mtdSlewingCorr"));
166  if(!mtdSlewingCorr) {
167  LOG_ERROR << "Unable to get the slewing correction parameters from data base" << endm;
168  return kStErr;
169  }
170  mtdSlewingCorr_st* slewingCorr = static_cast<mtdSlewingCorr_st*>(mtdSlewingCorr->GetTable());
171  for(int j=0; j<mNBackleg*mNModule*mNBinMax; j++)
172  {
173  Int_t backleg = j/100;
174  Int_t module = (j-backleg*100)/20;
175  Int_t bin = j-backleg*100-module*20;
176  mMtdTotEdge[backleg][module][bin] = slewingCorr->slewingBinEdge[j];
177  mMtdTotCorr[backleg][module][bin] = slewingCorr->slewingCorr[j];
178  }
179 
180  if(mUseTriggercut)
181  {
182  // extract trigger time window cuts
183  mDbDataSet = GetDataBase("Calibrations/mtd/mtdTriggerTimeCut");
184  St_mtdTriggerTimeCut* mtdTriggerTimeCut = static_cast<St_mtdTriggerTimeCut*>(mDbDataSet->Find("mtdTriggerTimeCut"));
185  if(!mtdTriggerTimeCut) {
186  LOG_ERROR << "Unable to get the trigger time window cuts from data base" << endm;
187  return kStErr;
188  }
189  mtdTriggerTimeCut_st* trigTimeCut = static_cast<mtdTriggerTimeCut_st*>(mtdTriggerTimeCut->GetTable());
190  for(int j=0; j<mNBackleg*mNModule; j++)
191  {
192  mTriggerLowEdge[j/5][j%5] = trigTimeCut->minTriggerTime[j];
193  mTriggerHighEdge[j/5][j%5] = trigTimeCut->maxTriggerTime[j];
194  }
195  }
196  LOG_INFO << "Done initializing calibration parameters from data base" << endm;
197  }
198  else
199  {
200  LOG_INFO << "Initializing calibration parameters from local files" << endm;
201  ifstream inData;
202 
203  //load T0 offset parameters from local file
204  if (mCalibFileT0.length()==0)
205  {
206  LOG_ERROR << "Please input the local file path for T0 offset parameters" << endm;
207  return kStErr;
208  }
209  if (mDebug) { LOG_INFO << " Local file for T0 offset : " << mCalibFileT0 << endm; }
210  inData.open(mCalibFileT0.c_str());
211  if(!inData.is_open())
212  {
213  LOG_ERROR << "Unable to get the T0 offset parameters from local file" <<endm;
214  LOG_ERROR << "Check if this file exists: " << mCalibFileT0.c_str() << endm;
215  return kStErr;
216  }
217  Int_t backlegId, moduleId, cellId;
218  Int_t nbin;
219  Double_t t0Corr;
220 
221  for(Int_t i=0;i<mNBackleg;i++)
222  {
223  for(Int_t j=0;j<mNModule;j++)
224  {
225  for(Int_t l=0;l<mNCell;l++)
226  {
227  inData>>backlegId>>moduleId>>cellId;
228  inData>>t0Corr;
229  if(backlegId>=1 && backlegId<=mNBackleg &&
230  moduleId>=1 && moduleId<=mNModule &&
231  cellId>=0 && cellId<=mNCell-1)
232  {
233  mMtdT0Corr[backlegId-1][moduleId-1][cellId]=t0Corr;
234  if (mDebug) { LOG_INFO << "mMtdT0Corr=" <<mMtdT0Corr[backlegId-1][moduleId-1][cellId] << endm;}
235  }
236  }
237  }
238  }
239  inData.close();
240 
241  //load dy offset parameters from local file
242  if (mCalibFileDy.length()==0)
243  {
244  LOG_ERROR << "Please input the local file path for Dy offset parameters" << endm;
245  return kStErr;
246  }
247  if (mDebug) { LOG_INFO << " Local file for Dy offset : " << mCalibFileDy << endm; }
248  inData.open(mCalibFileDy.c_str());
249  if(!inData.is_open())
250  {
251  LOG_ERROR << "Unable to get the Dy offset parameters from local file" <<endm;
252  LOG_ERROR << "Check if this file exists: " << mCalibFileDy.c_str() << endm;
253  return kStErr;
254  }
255  Double_t yCorr;
256 
257  for(Int_t i=0;i<mNBackleg;i++)
258  {
259  for(Int_t j=0;j<mNModule;j++)
260  {
261  inData>>backlegId>>moduleId;
262  inData>>yCorr;
263  if(backlegId>=1 && backlegId<=mNBackleg &&
264  moduleId>=1 && moduleId<=mNModule)
265  {
266  mMtdDyCorr[backlegId-1][moduleId-1]=yCorr;
267  if (mDebug) { LOG_INFO << "mMtdDyCorr=" <<mMtdDyCorr[backlegId-1][moduleId-1]<< endm;}
268  }
269  }
270  }
271  inData.close();
272 
273  //load dz offset parameters from local file
274  if (mCalibFileDz.length()==0)
275  {
276  LOG_ERROR << "Please input the local file path for Dz offset parameters" << endm;
277  return kStErr;
278  }
279  if (mDebug) { LOG_INFO << " Local file for Dz offset : " << mCalibFileDz << endm; }
280  inData.open(mCalibFileDz.c_str());
281  if(!inData.is_open())
282  {
283  LOG_ERROR << "Unable to get the Dz offset parameters from local file" <<endm;
284  LOG_ERROR << "Check if this file exists: " << mCalibFileDz.c_str() << endm;
285  return kStErr;
286  }
287  Double_t zCorr;
288 
289  for(Int_t i=0;i<mNBackleg;i++)
290  {
291  for(Int_t j=0;j<mNModule;j++)
292  {
293  for(Int_t l=0;l<mNCell;l++)
294  {
295  inData>>backlegId>>moduleId>>cellId;
296  inData>>zCorr;
297  if(backlegId>=1 && backlegId<=mNBackleg &&
298  moduleId>=1 && moduleId<=mNModule &&
299  cellId>=0 && cellId<=mNCell-1)
300  {
301  mMtdDzCorr[backlegId-1][moduleId-1][cellId]=zCorr;
302  if (mDebug) { LOG_INFO << "mMtdDzCorr=" <<mMtdDzCorr[backlegId-1][moduleId-1][cellId] << endm;}
303  }
304  }
305  }
306  }
307  inData.close();
308 
309  //load slewing correction parameters from local file
310  if (mCalibFileTot.length()==0)
311  {
312  LOG_ERROR << "Please input the local file path for slewing correction parameters" << endm;
313  return kStErr;
314  }
315  if (mDebug) { LOG_INFO << " Local file for slewing correction : " << mCalibFileTot << endm; }
316  inData.open(mCalibFileTot.c_str());
317  if(!inData.is_open())
318  {
319  LOG_ERROR << "Unable to get the slewing correction parameters from local file" <<endm;
320  LOG_ERROR << "Check if this file exists: " << mCalibFileTot.c_str() << endm;
321  return kStErr;
322  }
323 
324  for(Int_t i=0;i<mNBackleg;i++)
325  {
326  for(Int_t j=0;j<mNModule;j++)
327  {
328  inData>>backlegId>>moduleId;
329  inData>>nbin;
330  if (mDebug) { LOG_INFO << "BL,MOD=" <<backlegId<<","<<moduleId<<" ,bin="<<nbin << endm; }
331  if(nbin<0) continue;
332  for(Int_t k=0;k<nbin;k++)
333  {
334  if(backlegId>=1 && backlegId<=mNBackleg &&
335  moduleId>=1 && moduleId<=mNModule)
336  {
337  inData>>mMtdTotEdge[backlegId-1][moduleId-1][k];
338  if (mDebug) { LOG_INFO << "edge=" <<mMtdTotEdge[backlegId-1][moduleId-1][k] << endm; }
339  }
340  }
341 
342  for(Int_t k=0;k<nbin;k++)
343  {
344  if(backlegId>=1 && backlegId<=mNBackleg &&
345  moduleId>=1 && moduleId<=mNModule)
346  {
347  inData>>mMtdTotCorr[backlegId-1][moduleId-1][k];
348  if (mDebug) { LOG_INFO << "Corr=" <<mMtdTotCorr[backlegId-1][moduleId-1][k] << endm; }
349  }
350  }
351  }
352  }
353  inData.close();
354 
355  if(mUseTriggercut)
356  {
357  //load trigger time window cuts from local file
358  if (mCalibFileTrigger.length()==0)
359  {
360  LOG_ERROR << "Please input the local file path for trigger time window cuts" << endm;
361  return kStErr;
362  }
363  if (mDebug) { LOG_INFO << " Local file for trigger time window cuts : " << mCalibFileTrigger << endm; }
364 
365  inData.open(mCalibFileTrigger.c_str());
366  if(!inData.is_open())
367  {
368  LOG_ERROR << "Unable to get the trigger time window cuts from local file" <<endm;
369  LOG_ERROR << "Check if this file exists: " << mCalibFileTrigger.c_str() << endm;
370  return kStErr;
371  }
372  for(Int_t i=0;i<mNBackleg;i++)
373  {
374  for(Int_t j=0;j<mNModule;j++)
375  {
376  inData>>backlegId>>moduleId>>mTriggerLowEdge[i][j]>>mTriggerHighEdge[i][j];
377  if (mDebug) { LOG_INFO << "triglow = " <<mTriggerLowEdge[i][j]<<", trighigh = " <<mTriggerHighEdge[i][j] << endm; }
378  }
379  }
380  inData.close();
381  }
382  }
383 
384  return kStOK;
385 }
386 
387 
388 //_____________________________________________________________________________
389 Int_t StMtdCalibMaker::Make()
390 {
391  if (mDebug) { LOG_INFO << "StMtdCalibMaker::Maker: starting ..." << endm; }
392 
393  mStEvent = (StEvent*) GetInputDS("StEvent");
394  if(mStEvent)
395  {
396  if (mDebug) { LOG_INFO << "Running on StEvent ..." << endm; }
397  processStEvent();
398  return kStOK;
399  }
400  else
401  {
402  StMuDstMaker *muDstMaker = (StMuDstMaker*) GetMaker("MuDst");
403  if(muDstMaker)
404  {
405  mMuDst = muDstMaker->muDst();
406  if (mMuDst)
407  {
408  if (mDebug) { LOG_INFO << "Running on muDst ..." << endm; }
409  processMuDst();
410  return kStOK;
411  }
412  }
413  }
414 
415  LOG_WARN << "Neither StEvent nor muDst is available!" << endm;
416  return kStWarn;
417 }
418 
419 //_____________________________________________________________________________
420 void StMtdCalibMaker::processStEvent()
421 {
422  // Get collision start time from VPD
423  StBTofCollection *btofCollection = mStEvent->btofCollection();
424  if (!btofCollection) return;
425  StBTofHeader *bTofHeader = btofCollection->tofHeader();
426  if (!bTofHeader) return;
427  Double_t tStartTime = bTofHeader->tStart();
428  if (mDebug) { LOG_INFO << " Collision start time =" <<tStartTime << endm; }
429 
430  StMtdCollection *mtdCollection = mStEvent->mtdCollection();
431  if(!mtdCollection)
432  {
433  LOG_WARN << "No Mtd collection is available ..." << endm;
434  return;
435  }
436 
437  // Get trigger time recorded in THUB
438  StMtdHeader *mtdHeader = mtdCollection->mtdHeader();
439  Double_t mtdTriggerTime[2] = {0,0};
440  if (!mtdHeader)
441  {
442  LOG_WARN << "No mtd header ..." << endm;
443  return;
444  }
445  mtdTriggerTime[0] = 25.*(mtdHeader->triggerTime(0)&0xfff);
446  mtdTriggerTime[1] = 25.*(mtdHeader->triggerTime(1)&0xfff);
447 
448  StSPtrVecMtdHit& mtdHits = mtdCollection->mtdHits();
449  Int_t nMtdHits = mtdHits.size();
450 
451  // loop over all the MTD hits
452  for(Int_t i=0; i<nMtdHits; i++)
453  {
454  StMtdHit *aHit = mtdHits[i];
455  if(!aHit) continue;
456  if(aHit->idTruth()>0) continue; // Do not apply calibration on MC hits
457  Int_t backlegId = aHit->backleg();
458  Int_t moduleId = aHit->module();
459  Int_t cellId = aHit->cell();
460  pair<Double_t,Double_t> leadTime = aHit->leadingEdgeTime();//ns
461  if(backlegId<1 || backlegId>mNBackleg) continue;
462  if(moduleId<1 || moduleId>mNModule) continue;
463  if(cellId<0 || cellId>mNCell-1) continue;
464 
465  //trigger time cut, and fill trigger histo
466  Int_t thub = 2;
467  if(backlegId>=16 && backlegId<=30) thub = 1;
468  Double_t triggerTime = (leadTime.first+leadTime.second)/2-mtdTriggerTime[thub-1];
469  while (triggerTime<-51200) triggerTime+=51200;
470  if (mUseTriggercut &&
471  (triggerTime<mTriggerLowEdge[backlegId-1][moduleId-1]|| triggerTime>mTriggerHighEdge[backlegId-1][moduleId-1]))
472  continue;
473  if (mHisto) hTriggerTimeBL->Fill((backlegId-1)*5+(moduleId-1),triggerTime);
474 
475  //time of flight before calibration
476  Double_t timeOfFlight = (leadTime.first+leadTime.second)/2-tStartTime;
477  while (timeOfFlight>51200) timeOfFlight-=51200;
478  while (timeOfFlight<-51200) timeOfFlight+=51200;
479  if (mHisto) hTimeOfFlightModule->Fill((backlegId-1)*5+moduleId-1,timeOfFlight);
480 
481  // Get the calibration parameters for T0 and slewing
482  pair<Double_t,Double_t> totTime = aHit->tot();
483  Double_t tot = sqrt(totTime.first*totTime.second);//geometric mean of tot
484  Double_t AllCorr = mtdAllCorr(tot,backlegId,moduleId,cellId);
485  if (mDebug) { LOG_INFO << "TOT = "<< tot << endm; }
486 
487  // Apply the calibration parameters
488  timeOfFlight += AllCorr;
489  if (mHisto) hAllCorr->Fill(AllCorr);
490  if (mHisto) hTimeOfFlightCorr->Fill(timeOfFlight);
491  if (mHisto) hTimeOfFlightCorrModule->Fill((backlegId-1)*5+moduleId-1,timeOfFlight);
492 
493  // Save calibrated time-of-flight to mtd Pidtraits
494  // global track
495  StTrack *gTrack = aHit->associatedTrack();
496  if(!gTrack) continue;
497  StSPtrVecTrackPidTraits& traits = gTrack->pidTraits();
498  StMtdPidTraits *mtdPid = 0;
499  for(UInt_t it=0; it<traits.size(); it++)
500  {
501  if(traits[it]->detector()==kMtdId)
502  {
503  mtdPid = dynamic_cast<StMtdPidTraits*>(traits[it]);
504  break;
505  }
506  }
507  if(!mtdPid) continue;
508  mtdPid->setTimeOfFlight(timeOfFlight);
509  if (mDebug)
510  { LOG_INFO << "Time-of-flight = " <<timeOfFlight <<", Calibration parameter = "<<AllCorr << endm; }
511  //change dy in mtdPidTraits
512  Float_t dy = mtdPid->deltaY();
513  dy -= mMtdDyCorr[backlegId-1][moduleId-1];
514  if (mHisto) hDyModule->Fill((backlegId-1)*5+moduleId-1,dy);
515  mtdPid->setDeltaY(dy);
516  //change dz in mtdPidTraits
517  Float_t dz = mtdPid->deltaZ();
518  dz -= mMtdDzCorr[backlegId-1][moduleId-1][cellId];
519  mtdPid->setDeltaZ(dz);
520  if (mHisto) hDzModule->Fill((backlegId-1)*5+moduleId-1,dz);
521 
522 
523  // primary track
524  StTrack *pTrack = gTrack->node()->track(primary);
525  if(!pTrack) continue;
526  StMtdPidTraits *pmtdPid = 0;
527  StSPtrVecTrackPidTraits &ptraits = pTrack->pidTraits();
528  for(UInt_t it=0; it<ptraits.size(); it++)
529  {
530  if(ptraits[it]->detector()==kMtdId)
531  {
532  pmtdPid = dynamic_cast<StMtdPidTraits*>(ptraits[it]);
533  break;
534  }
535  }
536  if(!pmtdPid) continue;
537  pmtdPid->setTimeOfFlight(timeOfFlight);
538  pmtdPid->setDeltaY(dy);
539  pmtdPid->setDeltaZ(dz);
540 
541  }
542  return;
543 }
544 
545 
546 //_____________________________________________________________________________
547 void StMtdCalibMaker::processMuDst()
548 {
549  // Get primary vertex position
550  StMuEvent* mEvent = mMuDst->event();
551  if (!mEvent)
552  {
553  LOG_WARN << "No muEvent is available ..." << endm;
554  return;
555  }
556  Double_t tpcvz = mEvent->primaryVertexPosition().z();
557 
558  // Get collision start time & vertex posision from VPD
559  StBTofHeader *bTofHeader = mMuDst->btofHeader();
560  if(!bTofHeader)
561  {
562  LOG_WARN << "No TOF header is available ..." << endm;
563  return;
564  }
565  Double_t tStartTime = bTofHeader->tStart();
566  Double_t vpdvz = bTofHeader->vpdVz();
567  //if (tStartTime==-9999) return;
568  if (mHisto) hVertexzVsTpcz->Fill(vpdvz,tpcvz); //check vertex distribution
569 
570  // Get trigger time recorded in THUB
571  Int_t nhits = mMuDst->numberOfMTDHit();
572  if (mDebug) { LOG_INFO << "Total number of MTD hits is: " << nhits << endm; }
573  Double_t mtdTriggerTime[2] = {0,0};
574  StMuMtdHeader *muMtdHeader = mMuDst->mtdHeader();
575  if(!muMtdHeader)
576  {
577  LOG_WARN << "No MTD header is available ..." << endm;
578  return;
579  }
580  mtdTriggerTime[0] = 25.*(muMtdHeader->triggerTime(1)&0xfff);
581  mtdTriggerTime[1] = 25.*(muMtdHeader->triggerTime(2)&0xfff);
582  // loop over all the MTD hits
583  for(Int_t i=0;i<nhits;i++)
584  {
585  StMuMtdHit* aHit = mMuDst->mtdHit(i) ;
586  if(!aHit) continue;
587  if(aHit->idTruth()>0) continue; // Do not apply calibration on MC hits
588  Int_t backlegId = aHit->backleg();
589  Int_t moduleId = aHit->module();
590  Int_t cellId = aHit->cell();
591  pair<Double_t,Double_t> leadTime = aHit->leadingEdgeTime();//ns
592  pair<Double_t,Double_t> trailTime = aHit->trailingEdgeTime();
593  if(backlegId<1 || backlegId>mNBackleg) continue;
594  if(moduleId<1 || moduleId>mNModule) continue;
595  if(cellId<0 || cellId>mNCell-1) continue;
596 
597  //trigger time cut, and fill trigger histo
598  Int_t thub = 2;
599  if(backlegId>=16 && backlegId<=30) thub = 1;
600  Double_t triggerTime = (leadTime.first+leadTime.second)/2-mtdTriggerTime[thub-1];
601  while (triggerTime<-51200) triggerTime+=51200;
602  if (mUseTriggercut &&
603  (triggerTime<mTriggerLowEdge[backlegId-1][moduleId-1]|| triggerTime>mTriggerHighEdge[backlegId-1][moduleId-1]))
604  continue;
605  if (mHisto) hTriggerTimeBL->Fill((backlegId-1)*5+(moduleId-1),triggerTime);
606 
607  //time of flight before calibration
608  Double_t timeOfFlight = (leadTime.first+leadTime.second)/2-tStartTime;
609  while (timeOfFlight>51200) timeOfFlight-=51200;
610  while (timeOfFlight<-51200) timeOfFlight+=51200;
611  if (mHisto) hTimeOfFlightModule->Fill((backlegId-1)*5+moduleId-1,timeOfFlight);
612 
613  // Get the calibration parameters for T0 and slewing
614  Double_t totFirst = trailTime.first-leadTime.first;
615  if (totFirst<0) totFirst+=51200;
616  Double_t totSecond = trailTime.second-leadTime.second;
617  if (totSecond<0) totSecond+=51200;
618  Double_t tot = sqrt(totFirst*totSecond);
619  Double_t AllCorr = mtdAllCorr(tot,backlegId,moduleId,cellId);
620  if (mDebug) { LOG_INFO << "TOT = "<< tot << endm; }
621 
622  // Apply the calibration parameters
623  timeOfFlight += AllCorr;
624  if (mHisto) hAllCorr->Fill(AllCorr);
625  if (mHisto) hTimeOfFlightCorr->Fill(timeOfFlight);
626  if (mHisto) hTimeOfFlightCorrModule->Fill((backlegId-1)*5+moduleId-1,timeOfFlight);
627 
628 
629  // Save calibrated time-of-flight to mtd Pidtraits
630  // global track
631  Int_t index = aHit->index2Global();
632  if (mDebug) { LOG_INFO << "Associated global track index = "<< index << endm; }
633  if(index<0) continue;
634  StMuTrack *gTrack = mMuDst->globalTracks(index);
635  if(!gTrack) continue;
636  StMuMtdPidTraits pidMtd = gTrack->mtdPidTraits();
637  pidMtd.setTimeOfFlight (timeOfFlight);
638  if (mDebug) { LOG_INFO << "Time-of-flight = " <<timeOfFlight <<", Calibration parameter = "<<AllCorr << endm; }
639  //change dy in mtdPidTraits
640  Float_t dy = pidMtd.deltaY();
641  dy -= mMtdDyCorr[backlegId-1][moduleId-1];
642  pidMtd.setDeltaY(dy);
643  if (mHisto) hDyModule->Fill((backlegId-1)*5+moduleId-1,dy);
644  //change dz in mtdPidTraits
645  Float_t dz = pidMtd.deltaZ();
646  dz -= mMtdDzCorr[backlegId-1][moduleId-1][cellId];
647  pidMtd.setDeltaZ(dz);
648  if (mHisto) hDzModule->Fill((backlegId-1)*5+moduleId-1,dz);
649  gTrack->setMtdPidTraits(pidMtd);
650 
651  //primary track
652  index = aHit->index2Primary();
653  if(index<0) continue;
654  StMuTrack *pTrack = (StMuTrack *)mMuDst->array(muPrimary)->UncheckedAt(index);
655  if(!pTrack) continue;
656  StMuMtdPidTraits ppidMtd = pTrack->mtdPidTraits();
657  ppidMtd.setTimeOfFlight (timeOfFlight);
658  ppidMtd.setDeltaY(dy);
659  ppidMtd.setDeltaZ(dz);
660  pTrack->setMtdPidTraits(ppidMtd);
661  } // end mtd hits loop
662 
663  return;
664 }
665 
666 //_____________________________________________________________________________
667 Double_t StMtdCalibMaker::mtdAllCorr(const Double_t tot, const Int_t iBackleg, const Int_t iModule, const Int_t iCell)
668 {
669  // return the calibration parameters for T0 and slewing combined
670  Int_t backleg = iBackleg-1;
671  Int_t module = iModule-1;
672  Int_t cell = iCell;
673  Double_t mtdCorr = 0;
674  mtdCorr += mMtdT0Corr[backleg][module][cell]; // T0 offset
675 
676  // slewing correction
677  Int_t iTotBin = -1;
678  for(Int_t i=0;i<mNBinMax-1;i++)
679  {
680  if(tot>=mMtdTotEdge[backleg][module][i] && tot<mMtdTotEdge[backleg][module][i+1])
681  {
682  iTotBin = i;
683  break;
684  }
685  }
686  if(iTotBin>=0&&iTotBin<mNBinMax)
687  {
688  Double_t x1 = mMtdTotEdge[backleg][module][iTotBin];
689  Double_t x2 = mMtdTotEdge[backleg][module][iTotBin+1];
690  Double_t y1 = mMtdTotCorr[backleg][module][iTotBin];
691  Double_t y2 = mMtdTotCorr[backleg][module][iTotBin+1];
692  mtdCorr += y1 + (tot-x1)*(y2-y1)/(x2-x1);
693  }
694  else
695  {
696  if (mDebug) { LOG_INFO << " TOT out of range! totCorr set 0! Tot=" << tot <<" ,BL,MOD="<<iBackleg<<","<<iModule<< endm; }
697  }
698  return mtdCorr;
699 }
700 
701 //_____________________________________________________________________________
702 void StMtdCalibMaker::bookHistograms()
703 {
704  hTimeOfFlightCorr = new TH1D("hTimeOfFlightCorr","hTimeOfFlightCorr",4000,-1000,1000);
705  hAllCorr = new TH1D("hAllCorr","hAllCorr",1000,-500,500);
706  hTimeOfFlightModule =
707  new TH2F("hTimeOfFlightModule","hTimeOfFlightModule",mNBackleg*mNModule,0,mNBackleg*mNModule,1000,-500,500);
708  hTimeOfFlightCorrModule =
709  new TH2F("hTimeOfFlightCorrModule","hTimeOfFlightCorrModule",mNBackleg*mNModule,0,mNBackleg*mNModule,1000,-500,500);
710  hTriggerTimeBL = new TH2F("hTriggerTimeBL","hTriggerTimeBL",150,0,150,400,2700,3100);
711  hVertexzVsTpcz = new TH2F("hVertexzVsTpcz","hVertexzVsTpcz",800,-100,100,800,-100,100);
712  hTOFTimeOfFlightTray = new TH2F("hTOFTimeOfFlightTray","hTOFTimeOfFlightTray",120,0,120,1000,-1000,1000);
713  hDyModule = new TH2F("hDyModule","hDyModule",mNBackleg*mNModule,0,mNBackleg*mNModule,200,-50,50);
714  hDzModule = new TH2F("hDzModule","hDzModule",mNBackleg*mNModule,0,mNBackleg*mNModule,240,-60,60);
715  AddHist(hTimeOfFlightCorr);
716  AddHist(hAllCorr);
717  AddHist(hTimeOfFlightModule);
718  AddHist(hDyModule);
719  AddHist(hDzModule);
720  AddHist(hTimeOfFlightCorrModule);
721  AddHist(hTriggerTimeBL);
722  AddHist(hVertexzVsTpcz);
723  AddHist(hTOFTimeOfFlightTray);
724  return;
725 }
726 
727 
728 //
729 // $Id: StMtdCalibMaker.cxx,v 1.6 2016/07/28 14:17:44 marr Exp $
730 // $Log: StMtdCalibMaker.cxx,v $
731 // Revision 1.6 2016/07/28 14:17:44 marr
732 // Fix coverity check: check variable nbin
733 //
734 // Revision 1.5 2016/07/27 15:17:19 marr
735 // Fix coverity checks: check values of backleg, module, cell
736 //
737 // Revision 1.4 2016/07/27 14:21:15 marr
738 // Check not to apply calibration parameters to MC hits
739 //
740 // Revision 1.3 2015/02/04 22:35:24 marr
741 // Check the existance of the matched primary track
742 //
743 // Revision 1.2 2015/02/04 14:35:59 marr
744 // Added dy and dz correction for pidTraits
745 //
746 // Revision 1.1 2014/10/09 16:59:52 jeromel
747 // Reviewed version of the MtdCalibMaker - first commit
748 //
StMtdCalibMaker::Make
virtual Int_t Make()
Definition: StMtdCalibMaker.cxx:389
StMuDst::globalTracks
static TObjArray * globalTracks()
returns pointer to the global tracks list
Definition: StMuDst.h:303
StMuEvent::primaryVertexPosition
StThreeVectorF primaryVertexPosition(int vtx_id=-1) const
The StMuDst is supposed to be structured in &#39;physical events&#39;. Therefore there is only 1 primary vert...
Definition: StMuEvent.cxx:221
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
StMtdCalibMaker::mtdAllCorr
Double_t mtdAllCorr(const Double_t tot, const Int_t iBackleg, const Int_t iModule, const Int_t iCell)
Caculate all the correction value for time of flight.
Definition: StMtdCalibMaker.cxx:667
StMtdCalibMaker::Init
virtual Int_t Init()
Destructor.
Definition: StMtdCalibMaker.cxx:96
StMtdCalibMaker::~StMtdCalibMaker
virtual ~StMtdCalibMaker()
Default constructor.
Definition: StMtdCalibMaker.cxx:81
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
kStWarn
Definition: Stypes.h:42
kStOK
Definition: Stypes.h:40
StMuDst::btofHeader
static StBTofHeader * btofHeader()
returns pointer to the btofHeader - dongx
Definition: StMuDst.h:423
StMuDst::array
static TClonesArray * array(int type)
returns pointer to the n-th TClonesArray
Definition: StMuDst.h:262
kStErr
Definition: Stypes.h:44
TDataSet::Find
virtual TDataSet * Find(const char *path) const
Definition: TDataSet.cxx:362
Generated by   doxygen 1.8.5