00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066 #include <iostream>
00067 #include "StEvent.h"
00068 #include "StBTofCollection.h"
00069 #include "StBTofHit.h"
00070 #include "StBTofHeader.h"
00071 #include "StBTofPidTraits.h"
00072 #include "StEventTypes.h"
00073 #include "Stypes.h"
00074 #include "StThreeVectorD.hh"
00075 #include "StHelix.hh"
00076 #include "StTrackGeometry.h"
00077 #include "StTrackPidTraits.h"
00078 #include "StEventUtilities/StuRefMult.hh"
00079 #include "PhysicalConstants.h"
00080 #include "StPhysicalHelixD.hh"
00081 #include "tables/St_tofTOffset_Table.h"
00082 #include "tables/St_tofTotbCorr_Table.h"
00083 #include "tables/St_tofZbCorr_Table.h"
00084
00085 #include "tables/St_vertexSeed_Table.h"
00086
00087 #include "StBTofUtil/tofPathLength.hh"
00088 #include "StBTofUtil/StBTofHitCollection.h"
00089 #include "StBTofUtil/StBTofGeometry.h"
00090
00091 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00092 #include "StMuDSTMaker/COMMON/StMuDst.h"
00093 #include "StMuDSTMaker/COMMON/StMuBTofHit.h"
00094 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00095 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
00096 #include "StMuDSTMaker/COMMON/StMuBTofPidTraits.h"
00097
00098 #include "StBTofCalibMaker.h"
00099 #include "StVpdCalibMaker/StVpdCalibMaker.h"
00100
00101 #ifdef __ROOT__
00102 ClassImp(StBTofCalibMaker)
00103 #endif
00104
00106 const Double_t StBTofCalibMaker::VHRBIN2PS = 24.4140625;
00108 const Double_t StBTofCalibMaker::HRBIN2PS = 97.65625;
00110 const Double_t StBTofCalibMaker::TMAX = 51200.;
00112 const Double_t StBTofCalibMaker::VZDIFFCUT=6.;
00114 const Double_t StBTofCalibMaker::DCARCUT=1.;
00115 const Double_t StBTofCalibMaker::mC_Light = C_C_LIGHT/1.e9;
00116
00117
00118
00119 StBTofCalibMaker::StBTofCalibMaker(const char *name) : StMaker(name)
00120 {
00124 setVPDHitsCut(1,1);
00125 setOuterGeometry(true);
00126
00127 mEvent = 0;
00128 mBTofHeader = 0;
00129 mMuDst = 0;
00130 mZCalibType = NOTSET;
00131 mTotCalibType = NOTSET;
00132
00133 mSlewingCorr = kTRUE;
00134 mMuDstIn = kFALSE;
00135 mUseVpdStart = kTRUE;
00136
00137 setCreateHistoFlag(kFALSE);
00138 setHistoFileName("btofcalib.root");
00139
00140
00141 mInitFromFile = kFALSE;
00142
00143
00144 setCalibFileTot("/star/institutions/rice/calib/default/totCali_4DB.dat");
00145 setCalibFileZhit("/star/institutions/rice/calib/default/zCali_4DB.dat");
00146 setCalibFileT0("/star/institutions/rice/calib/default/t0_4DB.dat");
00147
00148
00149 StThreeVectorD MomFstPt(0.,0.,9999.);
00150 StThreeVectorD origin(0.,0.,0.);
00151 mBeamHelix = new StPhysicalHelixD(MomFstPt,origin,0.5*tesla,1.);
00152 }
00153
00154
00155 StBTofCalibMaker::~StBTofCalibMaker()
00156 { }
00157
00158
00159 void StBTofCalibMaker::resetPars()
00160 {
00161 memset(mTofTotEdge, 0, sizeof(mTofTotEdge));
00162 memset(mTofTotCorr, 0, sizeof(mTofTotCorr));
00163 memset(mTofZEdge, 0, sizeof(mTofZEdge) );
00164 memset(mTofZCorr, 0, sizeof(mTofZCorr) );
00165 memset(mTofTZero, 0, sizeof(mTofTZero) );
00166 }
00167
00168
00169 void StBTofCalibMaker::resetVpd()
00170 {
00171 memset(mVPDLeTime, 0, sizeof(mVPDLeTime));
00172 mTStart = -9999.;
00173 mTDiff = -9999.;
00174 mProjVtxZ = -9999.;
00175 mVPDVtxZ = -9999.;
00176 mVPDHitPatternEast = 0;
00177 mVPDHitPatternWest = 0;
00178 mNEast = 0;
00179 mNWest = 0;
00180 mValidStartTime = kFALSE;
00181 mNTzero = 0;
00182 }
00183
00184
00185 Int_t StBTofCalibMaker::Init()
00186 {
00187 resetPars();
00188 resetVpd();
00189
00190 mUseEventVertex = ! IAttr("UseProjectedVertex");
00191 if (mUseEventVertex) {
00192 LOG_INFO << "Use event vertex position." << endm;
00193 } else {
00194 LOG_INFO << "Use projected vertex position." << endm;
00195 }
00196
00197
00198 if(m_Mode) {
00199
00200 } else {
00201 setHistoFileName("");
00202 }
00203
00204 if (mHisto){
00205 bookHistograms();
00206 LOG_INFO << "Histograms are booked" << endm;
00207 if (mHistoFileName!="") {
00208 LOG_INFO << "Histograms will be stored in " << mHistoFileName.c_str() << endm;
00209 }
00210 }
00211
00212 return kStOK;
00213 }
00214
00215
00216 Int_t StBTofCalibMaker::InitRun(int runnumber)
00217 {
00218
00219
00221 Int_t val = initParameters(runnumber);
00222 if(val==kStOK) {
00223 mValidCalibPar = kTRUE;
00224 LOG_DEBUG << "Initialized valid calibration parameters." << endm;
00225 } else {
00226 mValidCalibPar = kFALSE;
00227 LOG_WARN << "No valid calibration parameters! " << endm;
00228 }
00229
00231 StVpdCalibMaker *vpdCalib = (StVpdCalibMaker *)GetMaker("vpdCalib");
00232 if(vpdCalib) {
00233 mUseVpdStart = vpdCalib->useVpdStart();
00234 if (mUseVpdStart) {LOG_INFO << "Found VPD Calibration Maker: use vpd for start timing" << endm;}
00235 else {LOG_INFO << "Found VPD Calibration Maker: vpd **NOT** used for start timing" << endm;}
00236 } else {
00237 mUseVpdStart = kFALSE;
00238 LOG_INFO << "NO VPD Calibration Maker found: vpd **NOT** used for start timing" << endm;
00239 }
00240
00242 if(!mUseVpdStart && !mUseEventVertex) {
00243 LOG_WARN << " Try to run calibration without VPD as the start time and DON'T use the event vertex! Wrong command! Exit!" << endm;
00244 return kStOK;
00245 }
00246
00247
00248 return kStOK;
00249
00250 }
00251
00252
00253 Int_t StBTofCalibMaker::initParameters(int runnumber)
00254 {
00257
00258 if (mInitFromFile){
00259 LOG_INFO << "Initializing calibration parameters from files" << endm;
00260 ifstream inData;
00261
00263 LOG_INFO << " - ToT : " << mCalibFileTot << endm;
00264 inData.open(mCalibFileTot.c_str());
00265 int trayId, moduleId, cellId, boardId;
00266 int nbin;
00267 int iCalibType;
00268 inData >> iCalibType;
00269
00270
00271
00272 mTotCalibType = calibtype(iCalibType);
00273
00274
00275 switch(mTotCalibType) {
00276
00277 case CELLCALIB:
00278 for(int i=0;i<mNTray;i++) {
00279 for(int j=0;j<mNModule;j++) {
00280 for(int l=0;l<mNCell;l++){
00281 inData>>trayId>>moduleId>>cellId;
00282 inData>>nbin;
00283 for(int k=0;k<=nbin;k++) inData>>mTofTotEdge[trayId-1][moduleId-1][cellId-1][k];
00284 for(int k=0;k<=nbin;k++) {
00285 inData>>mTofTotCorr[trayId-1][moduleId-1][cellId-1][k];
00286 if(k%10==0&&Debug()) {
00287 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;
00288 }
00289 }
00290 }
00291 }
00292 }
00293 break;
00294
00295 case MODULECALIB:
00296 for(int i=0;i<mNTray;i++) {
00297 for(int j=0;j<mNModule;j++) {
00298 inData>>trayId>>moduleId;
00299 inData>>nbin;
00300
00301 for(int k=0;k<=nbin;k++){
00302 inData>>mTofTotEdge[trayId-1][moduleId-1][0][k];
00303 for(int l=0;l<mNCell;l++){
00304 mTofTotEdge[trayId-1][moduleId-1][l][k]=mTofTotEdge[trayId-1][moduleId-1][0][k];
00305 }
00306 }
00307 for(int k=0;k<=nbin;k++) {
00308 inData>>mTofTotCorr[trayId-1][moduleId-1][0][k];
00309 for(int l=0;l<mNCell;l++){
00310 mTofTotCorr[trayId-1][moduleId-1][l][k]=mTofTotCorr[trayId-1][moduleId-1][0][k];
00311 }
00312
00313 if(k%10==0&&Debug()) {
00314 LOG_DEBUG << " ijlk= " << i << " " << j << " " << 0 << " " << k << " tot " << mTofTotEdge[trayId-1][moduleId-1][0][k] << " corr " << mTofTotCorr[trayId-1][moduleId-1][0][k] << endm;
00315 }
00316 }
00317
00318 }
00319 }
00320 break;
00321
00322 case BOARDCALIB:
00323 for(int i=0;i<mNTray;i++) {
00324 for(int j=0;j<mNTDIG;j++) {
00325 inData>>trayId>>boardId;
00326 inData>>nbin;
00327
00328 for(int k=0;k<=nbin;k++){
00329 inData>>mTofTotEdge[trayId-1][(boardId-1)*4][0][k];
00330 for(int m=0;m<4;m++){
00331 for(int l=0;l<mNCell;l++){
00332 mTofTotEdge[trayId-1][(boardId-1)*4+m][l][k] = mTofTotEdge[trayId-1][(boardId-1)*4][0][k];
00333 }
00334 }
00335 }
00336 for(int k=0;k<=nbin;k++) {
00337 inData>>mTofTotCorr[trayId-1][(boardId-1)*4][0][k];
00338 for(int m=0;m<4;m++){
00339 for(int l=0;l<mNCell;l++){
00340 mTofTotCorr[trayId-1][(boardId-1)*4+m][l][k] = mTofTotCorr[trayId-1][(boardId-1)*4][0][k];
00341 }
00342 }
00343 if(k%10==0&&Debug()) {
00344 LOG_DEBUG << " ijlk= " << i << " " << j*4 << " " << 0 << " " << k << " tot " << mTofTotEdge[trayId-1][(boardId-1)*4][0][k] << " corr " << mTofTotCorr[trayId-1][(boardId-1)*4][0][k] << endm;
00345 }
00346 }
00347 }
00348 }
00349 break;
00350
00351 default:
00352 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;
00353
00354 }
00355
00356 inData.close();
00357
00359 LOG_INFO << " - Zhit : " << mCalibFileZhit << endm;
00360 inData.open(mCalibFileZhit.c_str());
00361
00362 inData>>iCalibType;
00363 mZCalibType = calibtype(iCalibType);
00364
00365 switch(mZCalibType){
00366 case CELLCALIB:
00367 for(int i=0;i<mNTray;i++) {
00368 for(int j=0;j<mNModule;j++) {
00369 for(int l=0;l<mNCell;l++){
00370 inData>>trayId>>moduleId>>cellId;
00371 inData>>nbin;
00372 for(int k=0;k<=nbin;k++) inData>>mTofZEdge[trayId-1][moduleId-1][cellId-1][k];
00373 for(int k=0;k<=nbin;k++) {
00374 inData>>mTofZCorr[trayId-1][moduleId-1][cellId-1][k];
00375 if(k%10==0&&Debug()) {
00376 LOG_DEBUG << " ijlk= " << i << " " << j << " " << l << " " << k << " tot " << mTofZEdge[trayId-1][moduleId-1][cellId-1][k] << " corr " << mTofZCorr[trayId-1][moduleId-1][cellId-1][k] << endm;
00377 }
00378 }
00379 }
00380 }
00381 }
00382 break;
00383
00384 case MODULECALIB:
00385 for(int i=0;i<mNTray;i++) {
00386 for(int j=0;j<mNModule;j++) {
00387 inData>>trayId>>moduleId;
00388 inData>>nbin;
00389
00390 for(int k=0;k<=nbin;k++){
00391 inData>>mTofZEdge[trayId-1][moduleId-1][0][k];
00392 for(int l=0;l<mNCell;l++){
00393 mTofZEdge[trayId-1][moduleId-1][l][k]=mTofZEdge[trayId-1][moduleId-1][0][k];
00394 }
00395 }
00396 for(int k=0;k<=nbin;k++) {
00397 inData>>mTofZCorr[trayId-1][moduleId-1][0][k];
00398 for(int l=0;l<mNCell;l++){
00399 mTofZCorr[trayId-1][moduleId-1][l][k]=mTofZCorr[trayId-1][moduleId-1][0][k];
00400 }
00401 if(k%10==0&&Debug()) {
00402 LOG_DEBUG << " ijlk= " << i << " " << j << " " << 0 << " " << k << " tot " << mTofZEdge[trayId-1][moduleId-1][0][k] << " corr " << mTofZCorr[trayId-1][moduleId-1][0][k] << endm;
00403 }
00404 }
00405 }
00406 }
00407 break;
00408
00409 case BOARDCALIB:
00410 for(int i=0;i<mNTray;i++) {
00411 for(int j=0;j<mNTDIG;j++) {
00412 inData>>trayId>>boardId;
00413 inData>>nbin;
00414
00415 for(int k=0;k<=nbin;k++){
00416 inData>>mTofZEdge[trayId-1][(boardId-1)*4][0][k];
00417 for(int m=0;m<4;m++){
00418 for(int l=0;l<mNCell;l++){
00419 mTofZEdge[trayId-1][(boardId-1)*4+m][l][k] = mTofZEdge[trayId-1][(boardId-1)*4][0][k];
00420 }
00421 }
00422 }
00423 for(int k=0;k<=nbin;k++) {
00424 inData>>mTofZCorr[trayId-1][(boardId-1)*4][0][k];
00425 for(int m=0;m<4;m++){
00426 for(int l=0;l<mNCell;l++){
00427 mTofZCorr[trayId-1][(boardId-1)*4+m][l][k] = mTofZCorr[trayId-1][(boardId-1)*4][0][k];
00428 }
00429 }
00430 if(k%10==0&&Debug()) {
00431 LOG_DEBUG << " ijlk= " << i << " " << j*4 << " " << 0 << " " << k << " tot " << mTofTotEdge[trayId-1][(boardId-1)*4][0][k] << " corr " << mTofTotCorr[trayId-1][(boardId-1)*4][0][k] << endm;
00432 }
00433 }
00434 }
00435 }
00436
00437
00438 break;
00439 default:
00440 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;
00441
00442 }
00443 inData.close();
00444
00446 LOG_INFO << " - T0 : " << mCalibFileT0 << endm;
00447 inData.open(mCalibFileT0.c_str());
00448
00449 for(int i=0;i<mNTray;i++) {
00450 for(int j=0;j<mNModule;j++) {
00451 for(int k=0;k<mNCell;k++) {
00452 inData>>trayId>>moduleId>>cellId;
00453 inData>>mTofTZero[trayId-1][moduleId-1][cellId-1];
00454 int index = i*mNModule*mNCell+j*mNCell+k;
00455 if(index%100==0&&Debug()) {
00456 LOG_DEBUG << " ijk= " << i << " " << j << " " << k << " t0 " << mTofTZero[trayId-1][moduleId-1][cellId-1] << endm;
00457 }
00458 }
00459 }
00460 }
00461 inData.close();
00462 }
00463 else {
00464
00466 LOG_INFO << "Initializing calibration parameters from database" << endm;
00467
00468
00469 TDataSet *mDbDataSet = GetDataBase("Calibrations/tof/tofTotbCorr");
00470 St_tofTotbCorr* tofTotCorr = static_cast<St_tofTotbCorr*>(mDbDataSet->Find("tofTotbCorr"));
00471 if(!tofTotCorr) {
00472 LOG_ERROR << "unable to get tof TotbCorr table parameters" << endm;
00473
00474 return kStErr;
00475 }
00476 tofTotbCorr_st* totCorr = static_cast<tofTotbCorr_st*>(tofTotCorr->GetArray());
00477 Int_t numRows = tofTotCorr->GetNRows();
00478
00479 if(numRows!=mNTray*mNTDIG && numRows!=mNTray*mNModule*mNCell && numRows!=mNTray*mNModule) {
00480 LOG_WARN << " Mis-matched number of rows in tofTotbCorr table! " << numRows
00481 << " (exp:" << mNTray*mNTDIG << " or " << mNTray*mNModule*mNCell << " or "<< mNTray*mNModule << ")" << endm;
00482 }
00483
00484 LOG_INFO << " Number of rows read in: " << numRows << " for ToT correction" << endm;
00485
00486
00487 mTotCalibType = calibtype(numRows);
00488
00489 switch(mTotCalibType){
00490 case CELLCALIB:
00491 for (Int_t i=0;i<numRows;i++) {
00492 short trayId = totCorr[i].trayId;
00493 short moduleId = totCorr[i].moduleId;
00494 short cellId = totCorr[i].cellId;
00495
00496
00497 LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
00498 for(Int_t j=0;j<mNBinMax;j++) {
00499 if(trayId>0&&trayId<=mNTray&&moduleId>0&&moduleId<=mNModule&&cellId>0&&cellId<=mNCell){
00500 mTofTotEdge[trayId-1][moduleId-1][cellId-1][j] = totCorr[i].tot[j];
00501 mTofTotCorr[trayId-1][moduleId-1][cellId-1][j] = totCorr[i].corr[j];
00502 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; }
00503 }
00504 }
00505 }
00506 break;
00507
00508 case MODULECALIB:
00509 for (Int_t i=0;i<numRows;i++) {
00510 short trayId = totCorr[i].trayId;
00511 short moduleId = totCorr[i].moduleId;
00512 short cellId = totCorr[i].cellId;
00513
00514
00515 LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
00516 for(Int_t j=0;j<mNBinMax;j++) {
00517 if(trayId>0&&trayId<=mNTray&&moduleId>0&&moduleId<=mNModule){
00518 for(Int_t k=0;k<mNCell;k++){
00519 mTofTotEdge[trayId-1][moduleId-1][cellId-1+k][j] = totCorr[i].tot[j];
00520 mTofTotCorr[trayId-1][moduleId-1][cellId-1+k][j] = totCorr[i].corr[j];
00521 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; }
00522 }
00523 }
00524 }
00525 }
00526 break;
00527
00528 case BOARDCALIB:
00529 for (Int_t i=0;i<numRows;i++) {
00530 short trayId = totCorr[i].trayId;
00531 short moduleId = totCorr[i].moduleId;
00532 short cellId = totCorr[i].cellId;
00533 short boardId = (moduleId-1)/4+1;
00534
00535 LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
00536 for(Int_t j=0;j<mNBinMax;j++) {
00537 if(trayId>0&&trayId<=mNTray&&boardId>0&&boardId<=mNTDIG){
00538 for(Int_t k=0; k<4;k++){
00539 for(Int_t l=0;l<mNCell;l++){
00540 mTofTotEdge[trayId-1][moduleId-1+k][cellId-1+l][j] = totCorr[i].tot[j];
00541 mTofTotCorr[trayId-1][moduleId-1+k][cellId-1+l][j] = totCorr[i].corr[j];
00542 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; }
00543 }
00544 }
00545 }
00546 }
00547 }
00548 break;
00549
00550 default:
00551 LOG_WARN << "Number of rows in tofTotbCorr table mis-matched. "<<endl;
00552 }
00553
00554
00555
00556 mDbDataSet = GetDataBase("Calibrations/tof/tofZbCorr");
00557 St_tofZbCorr* tofZCorr = static_cast<St_tofZbCorr*>(mDbDataSet->Find("tofZbCorr"));
00558 if(!tofZCorr) {
00559 LOG_ERROR << "unable to get tof ZbCorr table parameters" << endm;
00560
00561 return kStErr;
00562 }
00563 tofZbCorr_st* zCorr = static_cast<tofZbCorr_st*>(tofZCorr->GetArray());
00564 numRows = tofZCorr->GetNRows();
00565
00566 if(numRows!=mNTray*mNTDIG && numRows!=mNTray*mNModule*mNCell && numRows !=mNTray*mNModule) {
00567 LOG_WARN << " Mis-matched number of rows in tofZbCorr table! " << numRows
00568 << " (exp:" << mNTray*mNTDIG << " or " << mNTray*mNModule*mNCell << " or "<< mNTray*mNModule << ")" << endm;
00569 }
00570 LOG_INFO << " Number of rows read in: " << numRows << " for Z correction" << endm;
00571
00572
00573 mZCalibType = calibtype(numRows);
00574
00575 switch(mZCalibType){
00576 case CELLCALIB:
00577 for (Int_t i=0;i<numRows;i++) {
00578 short trayId = totCorr[i].trayId;
00579 short moduleId = totCorr[i].moduleId;
00580 short cellId = totCorr[i].cellId;
00581
00582
00583 LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
00584 for(Int_t j=0;j<mNBinMax;j++) {
00585 if(trayId>0&&trayId<=mNTray&&moduleId>0&&moduleId<=mNModule&&cellId>0&&cellId<=mNCell) {
00586 mTofZEdge[trayId-1][moduleId-1][cellId-1][j] = zCorr[i].z[j];
00587 mTofZCorr[trayId-1][moduleId-1][cellId-1][j] = zCorr[i].corr[j];
00588 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; }
00589 }
00590 }
00591 }
00592 break;
00593
00594 case MODULECALIB:
00595 for (Int_t i=0;i<numRows;i++) {
00596 short trayId = totCorr[i].trayId;
00597 short moduleId = totCorr[i].moduleId;
00598 short cellId = totCorr[i].cellId;
00599 short boardId = (moduleId-1)/4+1;
00600
00601 LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
00602 for(Int_t j=0;j<mNBinMax;j++) {
00603 if(trayId>0&&trayId<=mNTray&&boardId>0&&moduleId<=mNModule) {
00604 for(Int_t k=0;k<mNCell;k++){
00605 mTofZEdge[trayId-1][moduleId-1][cellId-1+k][j] = zCorr[i].z[j];
00606 mTofZCorr[trayId-1][moduleId-1][cellId-1+k][j] = zCorr[i].corr[j];
00607 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; }
00608 }
00609 }
00610 }
00611 }
00612 break;
00613
00614 case BOARDCALIB:
00615 for (Int_t i=0;i<numRows;i++) {
00616 short trayId = totCorr[i].trayId;
00617 short moduleId = totCorr[i].moduleId;
00618 short cellId = totCorr[i].cellId;
00619 short boardId = (moduleId-1)/4+1;
00620
00621 LOG_DEBUG << " tray " << trayId << " module " << moduleId << " cell " << cellId << endm;
00622 for(Int_t j=0;j<mNBinMax;j++) {
00623 if(trayId>0&&trayId<=mNTray&&boardId>0&&boardId<=mNTDIG) {
00624 for(Int_t k=0;k<4;k++){
00625 for(Int_t l=0;l<mNCell;l++){
00626 mTofZEdge[trayId-1][moduleId-1+k][cellId-1+l][j] = zCorr[i].z[j];
00627 mTofZCorr[trayId-1][moduleId-1+k][cellId-1+l][j] = zCorr[i].corr[j];
00628 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; }
00629 }
00630 }
00631 }
00632 }
00633 }
00634 break;
00635
00636 default:
00637 LOG_WARN << "Number of rows in tofZbCorr table mis-matched. "<<endl;
00638 }
00639
00640 mDbDataSet = GetDataBase("Calibrations/tof/tofTOffset");
00641 St_tofTOffset* tofTOffset = static_cast<St_tofTOffset*>(mDbDataSet->Find("tofTOffset"));
00642 if(!tofTOffset) {
00643 LOG_ERROR << "unable to get tof TOffset table parameters" << endm;
00644
00645 return kStErr;
00646 }
00647 tofTOffset_st* tZero = static_cast<tofTOffset_st*>(tofTOffset->GetArray());
00648 numRows = tofTOffset->GetNRows();
00649
00650 LOG_DEBUG << " Number of rows read in: " << numRows << " for TOffset correction" << endm;
00651
00652 if(numRows!=mNTray) {
00653 LOG_WARN << " Mis-matched number of rows in tofTOffset table! " << numRows
00654 << " (exp:" << mNTray << ")" << endm;
00655
00656 }
00657 for (Int_t i=0;i<numRows;i++) {
00658 short trayId = tZero[i].trayId;
00659 LOG_DEBUG << " tray " << trayId << endm;
00660
00661 if(trayId>0&&trayId<=mNTray) {
00662 for(int j=0;j<mNTOF;j++) {
00663 mTofTZero[trayId-1][j/6][j%6] = tZero[i].T0[j];
00664 if(Debug()&&j%10==0) { LOG_DEBUG << " j=" << j << " T0 " << mTofTZero[trayId-1][j/6][j%6] << endm; }
00665 }
00666 }
00667 }
00668 }
00669
00670
00671
00672 double x0 = 0.;
00673 double y0 = 0.;
00674 double dxdz = 0.;
00675 double dydz = 0.;
00676
00677
00678 TDataSet* dbDataSet = this->GetDataBase("Calibrations/rhic/vertexSeed");
00679
00680 if (dbDataSet) {
00681 vertexSeed_st* vSeed = ((St_vertexSeed*) (dbDataSet->FindObject("vertexSeed")))->GetTable();
00682
00683 x0 = vSeed->x0;
00684 y0 = vSeed->y0;
00685 dxdz = vSeed->dxdz;
00686 dydz = vSeed->dydz;
00687 }
00688 else {
00689 LOG_WARN << "No database for beamline (Calibrations/rhic/vertexSeed)" << endm;
00690 }
00691
00692 LOG_INFO << "BeamLine Constraint: " << endm;
00693 LOG_INFO << "x(z) = " << x0 << " + " << dxdz << " * z" << endm;
00694 LOG_INFO << "y(z) = " << y0 << " + " << dydz << " * z" << endm;
00695
00696
00697
00698
00699
00700
00701 StThreeVectorD origin(x0,y0,0.0);
00702 double pt = 88889999;
00703 double nxy=::sqrt(dxdz*dxdz + dydz*dydz);
00704 if(nxy<1.e-5){
00705 LOG_WARN << "Beam line must be tilted!" << endm;
00706 nxy=dxdz=1.e-5;
00707 }
00708 double p0=pt/nxy;
00709 double px = p0*dxdz;
00710 double py = p0*dydz;
00711 double pz = p0;
00712 StThreeVectorD MomFstPt(px*GeV, py*GeV, pz*GeV);
00713 if(mBeamHelix) delete mBeamHelix;
00714 mBeamHelix = new StPhysicalHelixD(MomFstPt,origin,0.5*tesla,1.);
00715
00716
00717 return kStOK;
00718 }
00719
00720
00721 Int_t StBTofCalibMaker::FinishRun(int runnumber)
00722 {
00723 if(mBeamHelix) delete mBeamHelix;
00724 mBeamHelix = 0;
00725
00726 return kStOK;
00727 }
00728
00729
00730 Int_t StBTofCalibMaker::Finish()
00731 {
00732 if (mHistoFileName!="") writeHistograms();
00733 return kStOK;
00734 }
00735
00736
00737 Int_t StBTofCalibMaker::Make()
00738 {
00739 LOG_DEBUG << "StBTofCalibMaker::Maker: starting ..." << endm;
00740 if(!mValidCalibPar) {
00741 LOG_WARN << "No valid calibration parameters. Skip ... " << endm;
00742 return kStOK;
00743 }
00744
00745 initEvent();
00746 resetVpd();
00747 if(mUseVpdStart) {
00748 loadVpdData();
00749 }
00750
00751 if(!mMuDstIn) processStEvent();
00752 else processMuDst();
00753
00754 writeStartTime();
00755
00756 return kStOK;
00757 }
00758
00759
00760 void StBTofCalibMaker::processStEvent()
00761 {
00762
00763 if( !mEvent ) {LOG_WARN << "No StEvent" << endm; return;}
00764 if (!mEvent->btofCollection()) {LOG_WARN << "No BTOFCollection" << endm; return;}
00765 if (!mEvent->btofCollection()->hitsPresent()) {LOG_WARN << "No hits present" << endm; return;}
00766
00767 StBTofCollection *theTof = mEvent->btofCollection();
00768 StSPtrVecBTofHit &tofHits = theTof->tofHits();
00769 Int_t nhits = tofHits.size();
00770 LOG_INFO << " Fired TOF cells + upVPD tubes : " << nhits << endm;
00771
00772 if(mUseVpdStart) {
00773
00774 mEvtVtxZ = -9999.;
00775 mProjVtxZ = -9999.;
00776 float dcaRmin = 9999.;
00777
00778
00779 if(mUseEventVertex) {
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00795
00796 StPrimaryVertex *pVtx = mEvent->primaryVertex();
00797 if (pVtx){
00798 mEvtVtxZ = pVtx->position().z();
00799 }
00800 else {
00801 LOG_WARN << "No (default) primary vertex information for this (st-) event" << endm;
00802 };
00803
00804 tstart(mEvtVtxZ, &mTStart, &mTDiff);
00805
00806 } else {
00810 for(int i=0;i<nhits;i++) {
00811 StBTofHit *aHit = dynamic_cast<StBTofHit*>(tofHits[i]);
00812 if(!aHit) continue;
00813 int trayId = aHit->tray();
00814 if(trayId>0&&trayId<=mNTray) {
00815 StGlobalTrack *gTrack = dynamic_cast<StGlobalTrack*>(aHit->associatedTrack());
00816 if(!gTrack) continue;
00817 StTrackGeometry *theTrackGeometry = gTrack->geometry();
00818
00819 StThreeVectorD tofPos = theTrackGeometry->helix().at(theTrackGeometry->helix().pathLengths(*mBeamHelix).first);
00820 StThreeVectorD beamPos = mBeamHelix->at(theTrackGeometry->helix().pathLengths(*mBeamHelix).second);
00821 StThreeVectorD dcatof = tofPos - beamPos;
00822
00823 LOG_DEBUG<<" tofPos(x,y,z) = "<<tofPos.x()<<","<<tofPos.y()<<","<<tofPos.z()<<endm;
00824 LOG_DEBUG<<"beamPos(x,y,z) = "<<beamPos.x()<<","<<beamPos.y()<<","<<beamPos.z()<<endm;
00825 LOG_DEBUG<<" dca (x,y,z) = "<<dcatof.x()<<","<<dcatof.y()<<","<<dcatof.z()<<endm;
00826 LOG_DEBUG<<" 2D dca = "<<sqrt(pow(dcatof.x(),2)+pow(dcatof.y(),2))<<endm;
00827 LOG_DEBUG<<" 2D signed dca = "<<theTrackGeometry->helix().geometricSignedDistance(beamPos.x(),beamPos.y())<<endm;
00828
00830 if(fabs(tofPos.z()-mVPDVtxZ)>VZDIFFCUT) continue;
00831
00832 if(dcaRmin>dcatof.perp()) {
00833 mProjVtxZ = tofPos.z();
00834 dcaRmin = dcatof.perp();
00835 }
00836 }
00837 }
00838
00839 if(dcaRmin>DCARCUT) mProjVtxZ = -9999.;
00840 tstart(mProjVtxZ, &mTStart, &mTDiff);
00841
00842 }
00843
00844 } else {
00845
00846 StPrimaryVertex *pVtx = mEvent->primaryVertex();
00847 if(!pVtx) {
00848 LOG_WARN << " No primary vertex ... bye-bye" << endm;
00849 return;
00850 }
00851 mEvtVtxZ = pVtx->position().z();
00852
00853 tstart_NoVpd(theTof, pVtx, &mTStart);
00854
00855 }
00856
00857 LOG_INFO << "primVz = " << mEvtVtxZ << " projVz = " << mProjVtxZ << " vpdVz = " << mVPDVtxZ << endm;
00858 LOG_INFO << "Tstart = " << mTStart << " Tdiff = " << mTDiff << " NTzero = " << mNTzero << endm;
00859 LOG_INFO << "NWest = " << mNWest << " NEast = " << mNEast << " TdcSum West = " << mTSumWest << " East = " << mTSumEast << endm;
00860
00861
00862 if(mTStart<-1000.) {
00863 LOG_INFO << "No valid start time for this event. Skip ..." << endm;
00864 mValidStartTime = kFALSE;
00865 return;
00866 } else {
00867 mValidStartTime = kTRUE;
00868 }
00869
00870
00871
00872
00873
00874 for(int i=0;i<nhits;i++) {
00875 StBTofHit *aHit = dynamic_cast<StBTofHit*>(tofHits[i]);
00876 if(!aHit) continue;
00877 int trayId = aHit->tray();
00878 if(trayId<=0 || trayId>mNTray) continue;
00879
00880 StGlobalTrack *gTrack = dynamic_cast<StGlobalTrack*>(aHit->associatedTrack());
00881 if(!gTrack) {
00882 LOG_DEBUG << " No associated Track with this hit." << endm;
00883 continue;
00884 }
00885
00886 const StPtrVecTrackPidTraits& theTofPidTraits = gTrack->pidTraits(kTofId);
00887 if(!theTofPidTraits.size()) continue;
00888
00889 StTrackPidTraits *theSelectedTrait = theTofPidTraits[theTofPidTraits.size()-1];
00890 if(!theSelectedTrait) continue;
00891
00892 StBTofPidTraits *pidTof = dynamic_cast<StBTofPidTraits *>(theSelectedTrait);
00893 if(!pidTof) continue;
00894
00895 double tot = aHit->tot();
00896 double tdc = aHit->leadingEdgeTime();
00897 double tof = tdc - mTStart;
00898 Double_t zhit = pidTof->zLocal();
00899
00900 int moduleChan = (aHit->module()-1)*6 + (aHit->cell()-1);
00901 Double_t tofcorr = tofAllCorr(tof, tot, zhit, trayId, moduleChan);
00902 if(tofcorr<0.) {
00903 LOG_DEBUG << " Calibration failed! ... " << endm;
00904 continue;
00905 }
00906
00907 pidTof->setTimeOfFlight((Float_t)tofcorr);
00908
00910 StPrimaryTrack *pTrack = dynamic_cast<StPrimaryTrack *>(gTrack->node()->track(primary));
00911 StBTofPidTraits *ppidTof = 0;
00912 if(pTrack) {
00913 const StPtrVecTrackPidTraits& pTofPidTraits = pTrack->pidTraits(kTofId);
00914 if(!pTofPidTraits.size()) continue;
00915
00916 StTrackPidTraits *pSelectedTrait = pTofPidTraits[pTofPidTraits.size()-1];
00917 if(!pSelectedTrait) continue;
00918
00919 ppidTof = dynamic_cast<StBTofPidTraits *>(pSelectedTrait);
00920
00921 if(ppidTof && mUseEventVertex) {
00922 ppidTof->setTimeOfFlight((Float_t)tofcorr);
00923 }
00924 }
00925
00927 Double_t L = -9999.;
00928 Double_t ptot = -9999.;
00929 Bool_t doPID = kFALSE;
00930 if(mUseEventVertex) {
00931 if(!pTrack) {
00932 LOG_DEBUG << " The associated track is not a primary one. Skip PID calculation! " << endm;
00933 } else {
00934 StTrackGeometry *theTrackGeometry = pTrack->geometry();
00935 const StVertex *thisVertex = pTrack->vertex();
00936 if(!thisVertex) {
00937 LOG_DEBUG << " The associated track is not coming from any vertex. Skip PID calculation! " << endm;
00938 } else {
00939 StThreeVectorF primPos = thisVertex->position();
00940 L = tofPathLength(&primPos, &pidTof->position(), theTrackGeometry->helix().curvature());
00941 ptot = pTrack->geometry()->momentum().mag();
00942 doPID = kTRUE;
00943 }
00944 }
00945
00946 } else {
00947
00948 StTrackGeometry *theTrackGeometry = gTrack->geometry();
00949 StThreeVectorD tofPos = theTrackGeometry->helix().at(theTrackGeometry->helix().pathLengths(*mBeamHelix).first);
00950 StThreeVectorD dcatof = tofPos - mBeamHelix->at(theTrackGeometry->helix().pathLengths(*mBeamHelix).second);
00951 if(dcatof.perp()>DCARCUT) {
00952 LOG_DEBUG << " The projected position is far from beam line. Skip PID calculation! " << endm;
00953 } else if(fabs(tofPos.z()-mVPDVtxZ)>VZDIFFCUT) {
00954 LOG_DEBUG << " This track is not coming from the same VPD vertex! Skip PID calculation! " << endm;
00955 } else {
00956 L = tofPathLength(&tofPos, &pidTof->position(), theTrackGeometry->helix().curvature());
00957 ptot = gTrack->geometry()->momentum().mag();
00958 if(gTrack->dcaGeometry()) {
00959 ptot = gTrack->dcaGeometry()->momentum().mag();
00960 }
00961 doPID = kTRUE;
00962 }
00963
00964 }
00965
00966 if(!doPID) continue;
00967
00968 Double_t beta = L/(tofcorr*(C_C_LIGHT/1.e9));
00969
00970 Double_t b_e = ptot/sqrt(ptot*ptot+M_ELECTRON*M_ELECTRON);
00971 Double_t b_pi = ptot/sqrt(ptot*ptot+M_PION_PLUS*M_PION_PLUS);
00972 Double_t b_k = ptot/sqrt(ptot*ptot+M_KAON_PLUS*M_KAON_PLUS);
00973 Double_t b_p = ptot/sqrt(ptot*ptot+M_PROTON*M_PROTON);
00974
00975 float sigmae = -9999.;
00976 float sigmapi = -9999.;
00977 float sigmak = -9999.;
00978 float sigmap = -9999.;
00979 float res = 0.013;
00980 if(fabs(res)>1.e-5) {
00981 sigmae = (Float_t)((1./beta-1./b_e)/res);
00982 sigmapi = (Float_t)((1./beta-1./b_pi)/res);
00983 sigmak = (Float_t)((1./beta-1./b_k)/res);
00984 sigmap = (Float_t)((1./beta-1./b_p)/res);
00985 }
00986
00987 pidTof->setPathLength((Float_t)L);
00988 pidTof->setBeta((Float_t)beta);
00989 pidTof->setSigmaElectron(sigmae);
00990 pidTof->setSigmaPion(sigmapi);
00991 pidTof->setSigmaKaon(sigmak);
00992 pidTof->setSigmaProton(sigmap);
00993
00994 LOG_DEBUG << " storing BTofPidTraits for the global track" << endm;
00995
00996 if(mUseEventVertex) {
00997
00998 if(ppidTof) {
00999
01000 ppidTof->setPathLength((Float_t)L);
01001 ppidTof->setBeta((Float_t)beta);
01002 ppidTof->setSigmaElectron(sigmae);
01003 ppidTof->setSigmaPion(sigmapi);
01004 ppidTof->setSigmaKaon(sigmak);
01005 ppidTof->setSigmaProton(sigmap);
01006
01007 LOG_DEBUG << " storing BTofPidTraits for the primary track" << endm;
01008 }
01009 }
01010
01011 }
01012
01013 return;
01014 }
01015
01016
01017 void StBTofCalibMaker::processMuDst()
01018 {
01019 if(!mMuDst) {
01020 LOG_WARN << " No MuDst ... bye-bye" << endm;
01021 return;
01022 }
01023
01024 cleanCalibMuDst();
01025
01026 Int_t nhits = mMuDst->numberOfBTofHit();
01027 LOG_INFO << " Fired TOF cells + upVPD tubes : " << nhits << endm;
01028
01029 if(mUseVpdStart) {
01030
01031 mEvtVtxZ = -9999.;
01032 mProjVtxZ = -9999.;
01033 float dcaRmin = 9999.;
01034
01035
01036 if(mUseEventVertex) {
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01052
01053 StMuPrimaryVertex* pVtx = mMuDst->primaryVertex();
01054 if (pVtx){
01055 mEvtVtxZ = pVtx->position().z();
01056 }
01057 else {
01058 LOG_WARN << "No (default) primary vertex information for this (mudst) event" << endm;
01059 }
01060
01061 tstart(mEvtVtxZ, &mTStart, &mTDiff);
01062
01063 } else {
01067 for(int i=0;i<nhits;i++) {
01068 StMuBTofHit *aHit = (StMuBTofHit*)mMuDst->btofHit(i);
01069 if(!aHit) continue;
01070 int trayId = aHit->tray();
01071 if(trayId>0&&trayId<=mNTray) {
01072 StMuTrack *gTrack = aHit->globalTrack();
01073 if(!gTrack) continue;
01074
01075 StPhysicalHelixD thisHelix = gTrack->helix();
01076
01077 StThreeVectorD tofPos = thisHelix.at(thisHelix.pathLengths(*mBeamHelix).first);
01078 StThreeVectorD dcatof = tofPos - mBeamHelix->at(thisHelix.pathLengths(*mBeamHelix).second);
01079
01081 if(fabs(tofPos.z()-mVPDVtxZ)>VZDIFFCUT) continue;
01082
01083 if(dcaRmin>dcatof.perp()) {
01084 mProjVtxZ = tofPos.z();
01085 dcaRmin = dcatof.perp();
01086 }
01087 }
01088 }
01089
01090 if(dcaRmin>DCARCUT) mProjVtxZ = -9999.;
01091 tstart(mProjVtxZ, &mTStart, &mTDiff);
01092
01093 }
01094 } else {
01095
01096 StMuPrimaryVertex *pVtx = mMuDst->primaryVertex();
01097 if(!pVtx) {
01098 LOG_WARN << " No primary vertex ... bye-bye" << endm;
01099 return;
01100 }
01101 mEvtVtxZ = pVtx->position().z();
01102
01103 tstart_NoVpd(mMuDst, pVtx, &mTStart);
01104 }
01105
01106 LOG_INFO << "primVz = " << mEvtVtxZ << " projVz = " << mProjVtxZ << " vpdVz = " << mVPDVtxZ << endm;
01107 LOG_INFO << "Tstart = " << mTStart << " Tdiff = " << mTDiff << " NTzero = " << mNTzero << endm;
01108 LOG_INFO << "NWest = " << mNWest << " NEast = " << mNEast << " TdcSum West = " << mTSumWest << " East = " << mTSumEast << endm;
01109
01110
01111 if(mTStart<-1000.) {
01112 LOG_INFO << " No valid start time for this event. Skip ..." << endm;
01113 mValidStartTime = kFALSE;
01114 return;
01115 } else {
01116 mValidStartTime = kTRUE;
01117 }
01118
01119
01120
01121
01122
01123 for(int i=0;i<nhits;i++) {
01124 StMuBTofHit *aHit = (StMuBTofHit*)mMuDst->btofHit(i);
01125 if(!aHit) continue;
01126 int trayId = aHit->tray();
01127 if(trayId<=0 || trayId>mNTray) continue;
01128
01129 StMuTrack *gTrack = aHit->globalTrack();
01130 if(!gTrack) {
01131 LOG_DEBUG << " No associated Track with this hit." << endm;
01132 continue;
01133 }
01134
01135 StMuBTofPidTraits pidTof = gTrack->btofPidTraits();
01136
01137 double tot = aHit->tot();
01138 double tdc = aHit->leadingEdgeTime();
01139 while(tdc>TMAX) tdc -= TMAX;
01140 double tof = tdc - mTStart;
01141 Double_t zhit = pidTof.zLocal();
01142
01143 int moduleChan = (aHit->module()-1)*6 + (aHit->cell()-1);
01144 Double_t tofcorr = tofAllCorr(tof, tot, zhit, trayId, moduleChan);
01145 if(tofcorr<0.) {
01146 LOG_DEBUG << " Calibration failed! ... " << endm;
01147 continue;
01148 }
01149
01151 pidTof.setTimeOfFlight((Float_t)tofcorr);
01152
01154 StMuTrack *pTrack = aHit->primaryTrack();
01155 StMuBTofPidTraits ppidTof;
01156 if(pTrack) {
01157 ppidTof = pTrack->btofPidTraits();
01158 if(mUseEventVertex) ppidTof.setTimeOfFlight((Float_t)tofcorr);
01159 }
01160
01162 Double_t L = -9999.;
01163 Double_t ptot = -9999.;
01164 Bool_t doPID = kFALSE;
01165 if(mUseEventVertex) {
01166 if(!pTrack) {
01167 LOG_DEBUG << " The associated track is not a primary one. Skip PID calculation! " << endm;
01168 } else {
01169 int iv = pTrack->vertexIndex();
01170 StMuPrimaryVertex *thisVertex = mMuDst->primaryVertex(iv);
01171 if(!thisVertex) {
01172 LOG_DEBUG << " The associated track is not coming from any vertex. Skip PID calculation! " << endm;
01173 } else {
01174 StThreeVectorF primPos = thisVertex->position();
01175 StPhysicalHelixD thisHelix = pTrack->helix();
01176 L = tofPathLength(&primPos, &pidTof.position(), thisHelix.curvature());
01177 ptot = pTrack->momentum().mag();
01178 doPID = kTRUE;
01179 }
01180 }
01181
01182 } else {
01183
01184 StPhysicalHelixD gHelix = gTrack->helix();
01185 StThreeVectorD tofPos = gHelix.at(gHelix.pathLengths(*mBeamHelix).first);
01186 StThreeVectorD dcatof = tofPos - mBeamHelix->at(gHelix.pathLengths(*mBeamHelix).second);
01187 if(dcatof.perp()>DCARCUT) {
01188 LOG_DEBUG << " The projected position is far from beam line. Skip PID calculation! " << endm;
01189 } else if(fabs(tofPos.z()-mVPDVtxZ)>VZDIFFCUT) {
01190 LOG_DEBUG << " This track is not coming from the same VPD vertex! Skip PID calculation! " << endm;
01191 } else {
01192 L = tofPathLength(&tofPos, &pidTof.position(), gHelix.curvature());
01193 ptot = gTrack->momentum().mag();
01194 doPID = kTRUE;
01195 }
01196 }
01197
01198 if(doPID) {
01199 Double_t beta = L/(tofcorr*(C_C_LIGHT/1.e9));
01200
01201 Double_t b_e = ptot/sqrt(ptot*ptot+M_ELECTRON*M_ELECTRON);
01202 Double_t b_pi = ptot/sqrt(ptot*ptot+M_PION_PLUS*M_PION_PLUS);
01203 Double_t b_k = ptot/sqrt(ptot*ptot+M_KAON_PLUS*M_KAON_PLUS);
01204 Double_t b_p = ptot/sqrt(ptot*ptot+M_PROTON*M_PROTON);
01205
01206 float sigmae = -9999.;
01207 float sigmapi = -9999.;
01208 float sigmak = -9999.;
01209 float sigmap = -9999.;
01210 float res = 0.013;
01211 if(fabs(res)>1.e-5) {
01212 sigmae = (Float_t)((1./beta-1./b_e)/res);
01213 sigmapi = (Float_t)((1./beta-1./b_pi)/res);
01214 sigmak = (Float_t)((1./beta-1./b_k)/res);
01215 sigmap = (Float_t)((1./beta-1./b_p)/res);
01216 }
01217
01218 pidTof.setPathLength((Float_t)L);
01219 pidTof.setBeta((Float_t)beta);
01220 pidTof.setSigmaElectron(sigmae);
01221 pidTof.setSigmaPion(sigmapi);
01222 pidTof.setSigmaKaon(sigmak);
01223 pidTof.setSigmaProton(sigmap);
01224
01225 if(mUseEventVertex && pTrack) {
01226
01227 ppidTof.setPathLength((Float_t)L);
01228 ppidTof.setBeta((Float_t)beta);
01229 ppidTof.setSigmaElectron(sigmae);
01230 ppidTof.setSigmaPion(sigmapi);
01231 ppidTof.setSigmaKaon(sigmak);
01232 ppidTof.setSigmaProton(sigmap);
01233 }
01234 }
01235
01236 gTrack->setBTofPidTraits(pidTof);
01237 LOG_DEBUG << " storing BTofPidTraits for the global track" << endm;
01238
01239 if(mUseEventVertex && pTrack) {
01240 pTrack->setBTofPidTraits(ppidTof);
01241 LOG_DEBUG << " storing BTofPidTraits for the primary track" << endm;
01242 }
01243 }
01244
01245 return;
01246 }
01247
01248
01249 void StBTofCalibMaker::cleanCalibMuDst()
01250 {
01251 if(!mMuDst) return;
01252
01253 Int_t nPrimary = mMuDst->numberOfPrimaryTracks();
01254 Int_t nGlobal = mMuDst->numberOfGlobalTracks();
01255 for(int i=0;i<nPrimary;i++) {
01256 StMuTrack *pTrack = (StMuTrack *)mMuDst->primaryTracks(i);
01257 if(!pTrack) continue;
01258 StMuBTofPidTraits pid = pTrack->btofPidTraits();
01259 cleanCalib(pid);
01260 pTrack->setBTofPidTraits(pid);
01261 }
01262 for(int i=0;i<nGlobal;i++) {
01263 StMuTrack *gTrack = (StMuTrack *)mMuDst->globalTracks(i);
01264 if(!gTrack) continue;
01265 StMuBTofPidTraits pid = gTrack->btofPidTraits();
01266 cleanCalib(pid);
01267 gTrack->setBTofPidTraits(pid);
01268 }
01269 return;
01270 }
01271
01272
01273 void StBTofCalibMaker::cleanCalib(StMuBTofPidTraits& pid)
01274 {
01275 pid.setTimeOfFlight(-999.);
01276 pid.setPathLength(-999.);
01277 pid.setBeta(-999.);
01278 pid.setSigmaElectron(-999.);
01279 pid.setSigmaPion(-999.);
01280 pid.setSigmaKaon(-999.);
01281 pid.setSigmaProton(-999.);
01282 pid.setProbElectron(-999.);
01283 pid.setProbPion(-999.);
01284 pid.setProbKaon(-999.);
01285 pid.setProbProton(-999.);
01286 return;
01287 }
01288
01289
01290 void StBTofCalibMaker::initEvent()
01291 {
01292 if(mMuDstIn) {
01293 StMuDstMaker *mMuDstMaker = (StMuDstMaker *)GetMaker("MuDst");
01294 if(!mMuDstMaker) {
01295 LOG_WARN << " No MuDstMaker ... bye-bye" << endm;
01296 return;
01297 }
01298 mMuDst = mMuDstMaker->muDst();
01299 if(!mMuDst) {
01300 LOG_WARN << " No MuDst ... bye-bye" << endm;
01301 return;
01302 }
01303
01304 mBTofHeader = mMuDst->btofHeader();
01305 } else {
01306 mEvent = (StEvent *) GetInputDS("StEvent");
01307
01308 if(!mEvent) return;
01309 StBTofCollection *btofColl = mEvent->btofCollection();
01310 if(!btofColl) return;
01311 mBTofHeader = btofColl->tofHeader();
01312 }
01313
01314 return;
01315 }
01316
01317
01318 void StBTofCalibMaker::loadVpdData()
01319 {
01320 if(!mBTofHeader) return;
01321
01322 mTSumWest = 0;
01323 mTSumEast = 0;
01324 mVPDHitPatternWest = mBTofHeader->vpdHitPattern(west);
01325 mVPDHitPatternEast = mBTofHeader->vpdHitPattern(east);
01326 mNWest = mBTofHeader->numberOfVpdHits(west);
01327 mNEast = mBTofHeader->numberOfVpdHits(east);
01328 mVPDVtxZ = mBTofHeader->vpdVz();
01329
01330 for(int i=0;i<mNVPD;i++) {
01331 mVPDLeTime[i] = mBTofHeader->vpdTime(west, i+1);
01332 if(mVPDLeTime[i]>0.) mTSumWest += mVPDLeTime[i];
01333 if(Debug()) {
01334 LOG_DEBUG << " loading VPD West tubeId = " << i+1 << " time = " << mVPDLeTime[i] << endm;
01335 }
01336 }
01337
01338 for(int i=0;i<mNVPD;i++) {
01339 mVPDLeTime[i+mNVPD] = mBTofHeader->vpdTime(east, i+1);
01340 if(mVPDLeTime[i+mNVPD]>0.) mTSumEast += mVPDLeTime[i+mNVPD];
01341 if(Debug()) {
01342 LOG_DEBUG << " loading VPD East tubeId = " << i+1 << " time = " << mVPDLeTime[i+mNVPD] << endm;
01343 }
01344 }
01345
01346 return;
01347 }
01348
01349
01350 void StBTofCalibMaker::writeStartTime()
01351 {
01352 if(mBTofHeader) {
01353 mBTofHeader->setTStart(mTStart);
01354 mBTofHeader->setTDiff(mTDiff);
01355 mBTofHeader->setNTzero(mNTzero);
01356 }
01357
01358 return;
01359 }
01360
01361
01362 Double_t StBTofCalibMaker::tofAllCorr(const Double_t tof, const Double_t tot, const Double_t z, const Int_t iTray, const Int_t iModuleChan)
01363 {
01364 int tray = iTray;
01365 int module = iModuleChan/6 + 1;
01366 int cell = iModuleChan%6 + 1;
01367
01368 LOG_DEBUG << "\nStBTofCalibMaker::btofAllCorr: BTof calibrating...\n"
01369 << "\tDoing Calibration in BTOF Tray " << tray << " Module " << module << " Cell " << cell
01370 << "\n\tinput tof = " << tof
01371 << " TOT = " << tot << " Zlocal = " << z << endm;
01372
01373
01374 Double_t tofcorr = tof;
01375
01376 tofcorr -= mTofTZero[tray-1][module-1][cell-1];
01377
01378 LOG_DEBUG << "T0 correction: "<<mTofTZero[tray-1][module-1][cell-1]<<endm;
01379
01380 if(mSlewingCorr) {
01381 int iTotBin = -1;
01382 for(int i=0;i<mNBinMax-1;i++) {
01383 if(tot>=mTofTotEdge[tray-1][module-1][cell-1][i] && tot<mTofTotEdge[tray-1][module-1][cell-1][i+1]) {
01384 iTotBin = i;
01385 break;
01386 }
01387 }
01388 if(iTotBin>=0&&iTotBin<mNBinMax) {
01389 double x1 = mTofTotEdge[tray-1][module-1][cell-1][iTotBin];
01390 double x2 = mTofTotEdge[tray-1][module-1][cell-1][iTotBin+1];
01391 double y1 = mTofTotCorr[tray-1][module-1][cell-1][iTotBin];
01392 double y2 = mTofTotCorr[tray-1][module-1][cell-1][iTotBin+1];
01393 double dcorr = y1 + (tot-x1)*(y2-y1)/(x2-x1);
01394 LOG_DEBUG << "TOT correction: "<<dcorr<<endm;
01395
01396 tofcorr -= dcorr;
01397 } else {
01398 LOG_DEBUG << " TOT out of range! EXIT! " << endm;
01399 return -9999.;
01400 }
01401
01402 int iZBin = -1;
01403 for(int i=0;i<mNBinMax-1;i++) {
01404 if(z>=mTofZEdge[tray-1][module-1][cell-1][i] && z<mTofZEdge[tray-1][module-1][cell-1][i+1]) {
01405 iZBin = i;
01406 break;
01407 }
01408 }
01409 if(iZBin>=0&&iZBin<mNBinMax) {
01410 double x1 = mTofZEdge[tray-1][module-1][cell-1][iZBin];
01411 double x2 = mTofZEdge[tray-1][module-1][cell-1][iZBin+1];
01412 double y1 = mTofZCorr[tray-1][module-1][cell-1][iZBin];
01413 double y2 = mTofZCorr[tray-1][module-1][cell-1][iZBin+1];
01414 double dcorr = y1 + (z-x1)*(y2-y1)/(x2-x1);
01415
01416 tofcorr -= dcorr;
01417 LOG_DEBUG << "zHit correction: "<<dcorr<<endm;
01418
01419 } else {
01420 LOG_DEBUG << " Z out of range! EXIT! " << endm;
01421 return -9999.;
01422 }
01423
01424 }
01425
01426 LOG_DEBUG << " Corrected tof: tofcorr = " << tofcorr << endm;
01427 return tofcorr;
01428 }
01429
01430
01431 void StBTofCalibMaker::tstart(const Double_t vz, Double_t *tstart, Double_t *tdiff)
01432 {
01433 *tstart = -9999.;
01434 *tdiff = -9999.;
01435
01436 if(fabs(vz)>200.) {LOG_INFO << "tstart: vz too big" << endm; return;}
01437
01438 Double_t TSum = mTSumEast + mTSumWest;
01439
01440 if(mNEast+mNWest>0) {
01441 *tstart = (TSum-(mNEast-mNWest)*vz/(C_C_LIGHT/1.e9))/(mNEast+mNWest);
01442 }
01443 if ( mNEast>=mVPDEastHitsCut && mNWest>=mVPDWestHitsCut ) {
01444 *tdiff = (mTSumEast/mNEast - mTSumWest/mNWest)/2. - vz/(C_C_LIGHT/1.e9);
01445 }
01446
01447 return;
01448 }
01449
01450
01451 void StBTofCalibMaker::tstart_NoVpd(const StBTofCollection *btofColl, const StPrimaryVertex *pVtx, Double_t *tstart)
01452 {
01453 *tstart = -9999.;
01454 if(!btofColl) return;
01455
01456 const StSPtrVecBTofHit &tofHits = btofColl->tofHits();
01457 Int_t nCan = 0;
01458 Double_t tSum = 0.;
01459 Double_t t0[5000];
01460 memset(t0, 0., sizeof(t0));
01461 for(size_t i=0;i<tofHits.size();i++) {
01462 StBTofHit *aHit = dynamic_cast<StBTofHit*>(tofHits[i]);
01463 if(!aHit) continue;
01464 int trayId = aHit->tray();
01465 if(trayId>0&&trayId<=mNTray) {
01466 StGlobalTrack *gTrack = dynamic_cast<StGlobalTrack*>(aHit->associatedTrack());
01467 if(!gTrack) continue;
01468 StPrimaryTrack *pTrack = dynamic_cast<StPrimaryTrack*>(gTrack->node()->track(primary));
01469 if(!pTrack) continue;
01470 if(pTrack->vertex() != pVtx) continue;
01471 StThreeVectorF mom = pTrack->geometry()->momentum();
01472 double ptot = mom.mag();
01473
01474
01475 if(ptot<0.2 || ptot>0.6) continue;
01476
01477 static StTpcDedxPidAlgorithm PidAlgorithm;
01478 static StPionPlus* Pion = StPionPlus::instance();
01479 const StParticleDefinition* pd = pTrack->pidTraits(PidAlgorithm);
01480 double nSigPi = -999.;
01481 if(pd) {
01482 nSigPi = PidAlgorithm.numberOfSigma(Pion);
01483 }
01484
01485 if( fabs(nSigPi)>2.0 ) continue;
01486
01487 const StPtrVecTrackPidTraits& theTofPidTraits = pTrack->pidTraits(kTofId);
01488 if(!theTofPidTraits.size()) continue;
01489
01490 StTrackPidTraits *theSelectedTrait = theTofPidTraits[theTofPidTraits.size()-1];
01491 if(!theSelectedTrait) continue;
01492
01493 StBTofPidTraits *pidTof = dynamic_cast<StBTofPidTraits *>(theSelectedTrait);
01494 if(!pidTof) continue;
01495
01496 double tot = aHit->tot();
01497 double tof = aHit->leadingEdgeTime();
01498 double zhit = pidTof->zLocal();
01499
01500 int moduleChan = (aHit->module()-1)*6 + (aHit->cell()-1);
01501 Double_t tofcorr = tofAllCorr(tof, tot, zhit, trayId, moduleChan);
01502 if(tofcorr<0.) continue;
01503
01504 StThreeVectorF primPos = pVtx->position();
01505 StPhysicalHelixD helix = pTrack->geometry()->helix();
01506 double L = tofPathLength(&primPos, &pidTof->position(), helix.curvature());
01507 double tofPi = L*sqrt(M_PION_PLUS*M_PION_PLUS+ptot*ptot)/(ptot*(C_C_LIGHT/1.e9));
01508
01509 tSum += tofcorr - tofPi;
01510 t0[nCan] = tofcorr - tofPi;
01511 nCan++;
01512
01513 }
01514
01515 }
01516
01517 if(nCan<=0) {
01518 *tstart = -9999.;
01519 return;
01520 }
01521
01522 Int_t nTzero = nCan;
01523 if(nCan>1) {
01524 for(int i=0;i<nCan;i++) {
01525 double tdiff = t0[i] - (tSum-t0[i])/(nTzero-1);
01526 if(fabs(tdiff)>5.0) {
01527 tSum -= t0[i];
01528 nTzero--;
01529 }
01530 }
01531 }
01532
01533 mNTzero = nTzero;
01534
01535 *tstart = nTzero>0 ? tSum / nTzero : -9999.;
01536
01537 return;
01538 }
01539
01540
01541 void StBTofCalibMaker::tstart_NoVpd(const StMuDst *muDst, const StMuPrimaryVertex *pVtx, Double_t *tstart)
01542 {
01543 *tstart = -9999.;
01544 if(!muDst) return;
01545
01546 Int_t nBTofHits = muDst->numberOfBTofHit();
01547 Int_t nCan = 0;
01548 Double_t tSum = 0.;
01549 Double_t t0[5000];
01550 memset(t0, 0., sizeof(t0));
01551 for(int i=0;i<nBTofHits;i++) {
01552 StMuBTofHit *aHit = (StMuBTofHit*)muDst->btofHit(i);
01553 if(!aHit) continue;
01554 int trayId = aHit->tray();
01555 if(trayId>0&&trayId<=mNTray) {
01556 StMuTrack *gTrack = aHit->globalTrack();
01557 if(!gTrack) continue;
01558 StMuTrack *pTrack = aHit->primaryTrack();
01559 if(!pTrack) continue;
01560 StMuPrimaryVertex *aVtx = muDst->primaryVertex(pTrack->vertexIndex());
01561 if(aVtx != pVtx) continue;
01562 StThreeVectorF mom = pTrack->momentum();
01563 double ptot = mom.mag();
01564
01565
01566 if(ptot<0.2 || ptot>0.6) continue;
01567 double nSigPi = pTrack->nSigmaPion();
01568 if( fabs(nSigPi)>2. ) continue;
01569
01570 StMuBTofPidTraits pidTof = pTrack->btofPidTraits();
01571
01572 double tot = aHit->tot();
01573 double tof = aHit->leadingEdgeTime();
01574 double zhit = pidTof.zLocal();
01575
01576 int moduleChan = (aHit->module()-1)*6 + (aHit->cell()-1);
01577 Double_t tofcorr = tofAllCorr(tof, tot, zhit, trayId, moduleChan);
01578 if(tofcorr<0.) continue;
01579
01580 StThreeVectorF primPos = pVtx->position();
01581 StPhysicalHelixD helix = pTrack->helix();
01582 double L = tofPathLength(&primPos, &pidTof.position(), helix.curvature());
01583 double tofPi = L*sqrt(M_PION_PLUS*M_PION_PLUS+ptot*ptot)/(ptot*(C_C_LIGHT/1.e9));
01584
01585 tSum += tofcorr - tofPi;
01586 t0[nCan] = tofcorr - tofPi;
01587 nCan++;
01588
01589 }
01590
01591 }
01592
01593 if(nCan<=0) {
01594 *tstart = -9999.;
01595 return;
01596 }
01597
01598 Int_t nTzero = nCan;
01599 if(nCan>1) {
01600 for(int i=0;i<nCan;i++) {
01601 double tdiff = t0[i] - (tSum-t0[i])/(nTzero-1);
01602 if(fabs(tdiff)>5.0) {
01603 tSum -= t0[i];
01604 nTzero--;
01605 }
01606 }
01607 }
01608
01609 mNTzero = nTzero;
01610
01611 *tstart = nTzero>0 ? tSum / nTzero : -9999.;
01612
01613 return;
01614 }
01615
01616
01617 void StBTofCalibMaker::bookHistograms()
01618 {
01619 hEventCounter = new TH1D("eventCounter","eventCounter",20,0,20);
01620 }
01621
01622
01623 void StBTofCalibMaker::writeHistograms()
01624 {
01625
01626 TFile *theHistoFile = new TFile(mHistoFileName.c_str(), "RECREATE");
01627 LOG_INFO << "StBTofCalibMaker::writeHistograms()"
01628 << " histogram file " << mHistoFileName << endm;
01629
01630 theHistoFile->cd();
01631
01632 if(mHisto) {
01633 hEventCounter->Write();
01634 }
01635 return;
01636 }
