StRoot: StTofPool/StVpdAnalysisMaker/StVpdAnalysisMaker.cxx Source File

00001 /*******************************************************************
00002  *
00003  * $Id: StVpdAnalysisMaker.cxx,v 1.1 2008/09/02 18:27:21 dongx Exp $
00004  *
00005  * Author: Xin Dong
00006  *****************************************************************
00007  *
00008  * Description: Vpd analysis Maker to do the calibration for pVPD 
00009  *              hit pattern, VzVpd, Tdiff, Tstart
00010  *
00011  *****************************************************************
00012  *
00013  * $Log: StVpdAnalysisMaker.cxx,v $
00014  * Revision 1.1  2008/09/02 18:27:21  dongx
00015  * first release.
00016  * - Vpd analysis maker from MuDst to extract vz, Tstart, Tdiff etc.
00017  * - TPC primary vertex used for Tstart and Tdiff calculation
00018  *
00019  *
00020  *******************************************************************/
00021 #include <iostream>
00022 #include "TFile.h"
00023 #include "TProfile.h"
00024 #include "TNtuple.h"
00025 #include "TH1.h"
00026 #include "TH2.h"
00027 #include "TF1.h"
00028 
00029 #include "StEvent.h"
00030 #include "StEventTypes.h"
00031 #include "Stypes.h"
00032 #include "StMessMgr.h"
00033 #include "StThreeVectorD.hh"
00034 #include "phys_constants.h"
00035 #include "tables/St_tofTotCorr_Table.h"
00036 #include "tables/St_tofZCorr_Table.h"
00037 #include "tables/St_tofTOffset_Table.h"
00038 #include "tables/St_tofPhaseOffset_Table.h"
00039 
00040 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00041 #include "StMuDSTMaker/COMMON/StMuDst.h"
00042 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00043 
00044 #include "StVpdAnalysisMaker.h"
00045 
00046 #include "StMessMgr.h"
00047 #include "StMemoryInfo.hh"
00048 #include "StTimer.hh"
00049 
00050 #ifdef __ROOT__
00051 ClassImp(StVpdAnalysisMaker)
00052 #endif
00053 
00054 //_____________________________________________________________________________
00055 StVpdAnalysisMaker::StVpdAnalysisMaker(const char *name) : StMaker(name)
00056 {
00059   mTofINLCorr = 0;
00060   mDaqMap = 0;
00061   mVertexZ = -9999.;
00062 
00063   mEvent = 0;
00064   mMuDst = 0;
00065 
00066   setVPDHitsCut(1,1);
00067   mSlewingCorr = kTRUE;
00068   mValidCalibPar = kFALSE;
00069 
00070   resetCalibPars();
00071   resetVpdPars();
00072 }
00073 
00074 //_____________________________________________________________________________
00075 StVpdAnalysisMaker::~StVpdAnalysisMaker()
00076 { /* noop */ }
00077 
00078 //_____________________________________________________________________________
00079 void StVpdAnalysisMaker::resetCalibPars()
00080 {
00081   for(int i=0;i<2*mNVPD;i++) {
00082     for(int j=0;j<mNBinMax;j++) {
00083       mVPDTotEdge[i][j] = 0.0;
00084       mVPDTotCorr[i][j] = 0.0;
00085     }
00086     mVPDTZero[i] = 0.0;
00087     mVPDLeTime[i] = 0.0;
00088     mVPDTot[i] = 0.0;
00089   }
00090   mPhaseOffset8 = 0.;
00091 }
00092 
00093 //_____________________________________________________________________________
00094 void StVpdAnalysisMaker::resetVpdPars()
00095 {
00096   mTSumEast = 0.;
00097   mTSumWest = 0.;
00098   mVPDHitPatternEast = 0;
00099   mVPDHitPatternWest = 0;
00100   mNEast = 0;
00101   mNWest = 0;
00102   mVPDVtxZ = -9999.;
00103   mTDiff   = -9999.;
00104   mTStart  = -9999.;
00105   mVertexZ = -9999.;
00106 }
00107 
00108 //____________________________________________________________________________
00109 Int_t StVpdAnalysisMaker::Init()
00110 {
00111   return kStOK;
00112 }
00113 
00114 //____________________________________________________________________________
00115 Int_t StVpdAnalysisMaker::InitRun(int runnumber)
00116 {
00117   // tof run configurations
00118   mYear8 = (runnumber>9000000&&runnumber<10000000);
00119 
00120   mTofINLCorr = new StTofINLCorr();
00121   mDaqMap = new StTofrDaqMap();
00122   if(mYear8) {
00123     mTofINLCorr->setNValidTrays(mNValidTrays_Run8);
00124     mTofINLCorr->initFromDbase(this);
00125     gMessMgr->Info("","OS") << " Initialize INL table for run 8 ... " << endm;
00126 
00127     mDaqMap->setNValidTrays(mNValidTrays_Run8);
00128     mDaqMap->initFromDbaseGeneral(this);
00129     gMessMgr->Info("","OS") << " Initialize Daq map for run 8 ... " << endm;
00130   }
00131 
00132 
00133   Int_t val = kStOK;
00134   val = initParameters(runnumber);
00135   if(val==kStOK) {
00136     mValidCalibPar = kTRUE;
00137   } else {
00138     mValidCalibPar = kFALSE;
00139   }
00140 
00141   if(mValidCalibPar) {
00142     gMessMgr->Info(" ==> Good! Valid cali parameters! ","OS");
00143   } else {
00144     gMessMgr->Info(" ==> No valid cali parameters! ","OS");
00145   }
00146 
00147   return kStOK;
00148 
00149 }
00150 
00151 //_____________________________________________________________________________
00152 Int_t StVpdAnalysisMaker::initParameters(int runnumber)
00153 {
00156   LOG_INFO << "   -- retrieving run parameters from Calibrations_tof" << endm;
00157   TDataSet *mDbDataSet = GetDataBase("Calibrations/tof");
00158   if (!mDbDataSet){
00159     gMessMgr->Error("unable to get TOF run parameters","OS");
00160     return kStErr;
00161   }
00162   
00163   if(mYear8) {
00164     gMessMgr->Info("","OS") << "     loading parameters for Run VIII" << endm;
00165 
00166     // read tofTotCorr table
00167     St_tofTotCorr* tofTotCorr = static_cast<St_tofTotCorr*>(mDbDataSet->Find("tofTotCorr"));
00168     if(!tofTotCorr) {
00169       gMessMgr->Error("unable to get tof TotCorr table parameters","OS");
00170       //    assert(tofTotCorr);
00171       return kStErr;
00172     }
00173     tofTotCorr_st* totCorr = static_cast<tofTotCorr_st*>(tofTotCorr->GetArray());
00174     Int_t numRows = tofTotCorr->GetNRows();
00175 
00176     if(numRows!=mNTray8*mNTDIG+mNVPD*2) {
00177       gMessMgr->Warning("","OS") << " Mis-matched number of rows in tofTotCorr table! " << endm;
00178       //      return kStErr;
00179     }
00180 
00181     if(Debug()) gMessMgr->Info("","OS") << " Number of rows read in: " << numRows << " for ToT correction" << endm;
00182 
00183     for (Int_t i=0;i<mNTray8*mNTDIG+mNVPD*2;i++) {
00184       short trayId = totCorr[i].trayId;
00185       short moduleId = totCorr[i].moduleId;
00186       short boardId = (moduleId-1)/4+1;      // used for trays
00187       short cellId = totCorr[i].cellId;      // used for vpds
00188 
00189       int index = (trayId-120-1)*mNVPD+(cellId-1);   // used for vpd index
00190 
00191       if(Debug()) gMessMgr->Info("","OS") << " tray " << trayId << " board " << boardId << " cell " << cellId << endm;
00192       if(trayId!=mWestVpdTrayId&&trayId!=mEastVpdTrayId) continue;
00193       for(Int_t j=0;j<mNBinMax;j++) {
00194         if(trayId==mWestVpdTrayId||trayId==mEastVpdTrayId) { // upVPD east west
00195           mVPDTotEdge[index][j] = totCorr[i].tot[j];
00196           mVPDTotCorr[index][j] = totCorr[i].corr[j];
00197         } 
00198       } // end j 0->mNBinMax
00199     } // end i 0->numRows
00200 
00201     // read tofPhaseOffset table
00202     St_tofPhaseOffset* tofPhaseOffset = static_cast<St_tofPhaseOffset*>(mDbDataSet->Find("tofPhaseOffset"));
00203     if(!tofPhaseOffset) {
00204       gMessMgr->Error("unable to get tof PhaseOffset table parameters","OS");
00205       //    assert(tofPhaseOffset);
00206       return kStErr;
00207     }
00208     tofPhaseOffset_st* tPhaseDiff = static_cast<tofPhaseOffset_st*>(tofPhaseOffset->GetArray());
00209 
00210     mPhaseOffset8 = tPhaseDiff[0].T0[0] + tPhaseDiff[0].T0[1]; // run 8, only e/w difference, to double precision
00211     if(Debug()) gMessMgr->Info("","OS") << " PhaseOffset = " << mPhaseOffset8 << endm;
00212   }
00213 
00214   return kStOK;
00215 }
00216 //____________________________________________________________________________
00217 Int_t StVpdAnalysisMaker::FinishRun(int runnumber)
00218 {
00219   LOG_INFO << "StVpdAnalysisMaker -- cleaning up (FinishRun)" << endm;
00220 
00221   if(mDaqMap) delete mDaqMap;
00222   mDaqMap = 0;
00223 
00224   if(mTofINLCorr) delete mTofINLCorr;
00225   mTofINLCorr = 0;
00226   
00227   return kStOK;
00228 }
00229 
00230 //_____________________________________________________________________________
00231 Int_t StVpdAnalysisMaker::Finish()
00232 {
00233   return kStOK;
00234 }
00235 
00236 //_____________________________________________________________________________
00237 Int_t StVpdAnalysisMaker::Make()
00238 {
00239   gMessMgr->Info(" StVpdAnalysisMaker::Maker: starting ...","OS");
00240 
00241   resetVpdPars();
00242 
00243   Int_t iret = kStOK;
00244   if(mYear8) {
00245     iret = processEventYear8();
00246   }
00247   return kStOK;
00248 }
00249 
00250 //____________________________________________________________________________
00251 Int_t StVpdAnalysisMaker::processEventYear8()
00252 {
00253   mMuDstMaker = (StMuDstMaker *)GetMaker("MuDst");
00254   if(!mMuDstMaker) {
00255     LOG_INFO << " StVpdAnalysisMaker -- No MuDstMaker ... bye-bye" << endm;
00256     return kStOK;
00257   }
00258 
00259   mMuDst = mMuDstMaker->muDst();
00260   if(!mMuDst) {
00261     LOG_INFO << " StVpdAnalysisMaker -- No MuDst ... bye-bye" << endm;
00262     return kStOK;
00263   }
00264 
00265   mMuEvent = mMuDst->event();
00266   if(mMuEvent) {
00267     mVertexZ = mMuEvent->primaryVertexPosition().z();  
00268     if( fabs(mMuEvent->primaryVertexPosition().x())<1.e-4 &&
00269         fabs(mMuEvent->primaryVertexPosition().y())<1.e-4 &&
00270         fabs(mMuEvent->primaryVertexPosition().z())<1.e-4 )
00271     mVertexZ = -9999;
00272   }
00276   mEvent = new StEvent();
00277   mTofCollection = new StTofCollection();
00278   mEvent->setTofCollection(mTofCollection);
00279   int nTofRawData = mMuDst->numberOfTofRawData();
00280   for(int i=0;i<nTofRawData;i++) {
00281     StTofRawData *aRawData;
00282     StTofRawData *aMuRawData = (StTofRawData *)mMuDst->tofRawData(i);
00283     if(aMuRawData) {
00284       unsigned short tray = aMuRawData->tray();
00285       if(tray!=mWestVpdTrayId&&tray!=mEastVpdTrayId) continue;
00286       unsigned short leteFlag = aMuRawData->leteFlag();
00287       unsigned short channel = aMuRawData->channel();  
00288       unsigned int tdc = aMuRawData->tdc();
00289       unsigned int triggertime = aMuRawData->triggertime();
00290       unsigned short quality = aMuRawData->quality();
00291       aRawData = new StTofRawData(leteFlag,tray,channel,tdc,triggertime,quality);
00292     }
00293     mTofCollection->addRawData(aRawData);
00294   }
00295 
00296   mSortTofRawData = new StSortTofRawData(mTofCollection, mDaqMap);
00297 
00298   //-------------------------------------------------
00299   // upVPD calibration and calculate the start timing
00300   //-------------------------------------------------
00301   for(int i=0;i<2*mNVPD;i++) {
00302     mVPDLeTime[i] = -999.;
00303     mVPDTot[i] = -999.;
00304   }
00305   for(int ivpd=0;ivpd<2;ivpd++) { // west and east sides
00306     int trayId = (ivpd==0) ? mWestVpdTrayId : mEastVpdTrayId;
00307     int ewId = trayId - 120; 
00308 
00309     IntVec validtube = mSortTofRawData->GetValidChannel(trayId);
00310     if(Debug()) gMessMgr->Info("","OS") << " Number of fired hits on tray(vpd) " << trayId << " = " << validtube.size() << endm;
00311 
00312     if(!validtube.size()) continue;
00313     for(int i=0;i<mNVPD;i++) {
00314       int tubeId = i+1;
00315       int lechan = (ewId==1) ? mDaqMap->WestPMT2TDIGLeChan(tubeId) : mDaqMap->EastPMT2TDIGLeChan(tubeId);
00316       int techan = (ewId==1) ? mDaqMap->WestPMT2TDIGTeChan(tubeId) : mDaqMap->EastPMT2TDIGTeChan(tubeId);
00317       IntVec leTdc = mSortTofRawData->GetLeadingTdc(trayId, lechan, kTRUE);
00318       IntVec teTdc = mSortTofRawData->GetTrailingTdc(trayId, lechan, kTRUE);
00319 
00320       if(!leTdc.size() || !teTdc.size()) continue;
00321 
00322       int letdc = leTdc[0];
00323       int bin = (int)letdc&0x3ff;
00324       double letdc_f = letdc + mTofINLCorr->getVpdINLCorr(ewId, lechan, bin);
00325       double letime = letdc_f*VHRBIN2PS / 1000.;
00326 
00327       int tetdc = teTdc[0];
00328       bin = (int)tetdc&0x3ff;
00329       double tetdc_f = tetdc + mTofINLCorr->getVpdINLCorr(ewId, techan, bin);
00330       double tetime = tetdc_f*VHRBIN2PS / 1000.;
00331 
00332       mVPDLeTime[(trayId-121)*mNVPD+(tubeId-1)] = letime;
00333       mVPDTot[(trayId-121)*mNVPD+(tubeId-1)] = tetime - letime;
00334     }
00335   }
00336   tsum(mVPDTot, mVPDLeTime);
00337 
00338   gMessMgr->Info("","OS") << " NWest = " << mNWest << " NEast = " << mNEast << " TdcSum West = " << mTSumWest << " East = " << mTSumEast << endm;
00339 
00341   mTofCollection->setVpdEast(mVPDHitPatternEast);
00342   mTofCollection->setVpdWest(mVPDHitPatternWest);
00343   mTofCollection->setTdiff(mTDiff);
00344   mTofCollection->setVzVpd(mVPDVtxZ);
00345   mTofCollection->setTstart(mTStart);
00346 
00347   mNWest = mTofCollection->numberOfVpdWest();
00348   mNEast = mTofCollection->numberOfVpdEast();
00349 
00350   gMessMgr->Info("","OS") << " TofCollection: NWest = " << mVPDHitPatternWest << " NEast = " << mVPDHitPatternEast << endm;
00351   gMessMgr->Info("","OS") << "Tdiff = " << mTDiff <<" vpd vz = " << mVPDVtxZ << " tstart = " << mTStart << endm;
00352  
00353   return kStOK;
00354 }
00355 
00356 //_____________________________________________________________________________
00357 void StVpdAnalysisMaker::tsum(const Double_t *tot, const Double_t *time)
00358 {
00360   mTSumEast = 0.;
00361   mTSumWest = 0.;
00362   mNEast = 0;
00363   mNWest = 0;
00364   mVPDHitPatternWest = 0;
00365   mVPDHitPatternEast = 0;
00366 
00367   // West
00368   for(int i=0;i<mNVPD;i++) {
00369     if( time[i]>0. && tot[i]>0. ) {
00370       // west no offset
00371       if(mSlewingCorr) {   // a switch
00372         int ibin = -1;
00373         for(int j=0;j<mNBinMax-1;j++) {
00374           if(tot[i]>=mVPDTotEdge[i][j] && tot[i]<mVPDTotEdge[i][j+1]) {
00375             ibin = j;
00376             break;
00377           }
00378         }
00379         if(ibin>=0&&ibin<mNBinMax) {
00380           mNWest++;
00381 
00382           double x1 = mVPDTotEdge[i][ibin];
00383           double x2 = mVPDTotEdge[i][ibin+1];
00384           double y1 = mVPDTotCorr[i][ibin];
00385           double y2 = mVPDTotCorr[i][ibin+1];
00386           double dcorr = y1 + (tot[i]-x1)*(y2-y1)/(x2-x1);
00387 
00388           mVPDLeTime[i] = time[i] - dcorr;
00389           mTSumWest += mVPDLeTime[i];
00390 
00391           mVPDHitPatternWest |= 1<<i;
00392         }
00393       } else {
00394         gMessMgr->Info("","OS") << " Vpd West tube " << i+1 << " TOT out of range!" << endm;
00395         // out of range, remove this hit
00396       }
00397     }
00398       
00399     // East
00400     if( time[i+mNVPD]>0. && tot[i+mNVPD]>0. ) {
00401       double tmp = time[i+mNVPD] - mPhaseOffset8;
00402       while(tmp>TMAX) tmp -= TMAX;
00403       if(mSlewingCorr) {   // a switch
00404         int ibin = -1;
00405         for(int j=0;j<mNBinMax-1;j++) {
00406           if(tot[i+mNVPD]>=mVPDTotEdge[i+mNVPD][j] && tot[i+mNVPD]<mVPDTotEdge[i+mNVPD][j+1]) {
00407             ibin = j;
00408             break;
00409           }
00410         }
00411         if(ibin>=0&&ibin<mNBinMax) {
00412           mNEast++;
00413 
00414           double x1 = mVPDTotEdge[i+mNVPD][ibin];
00415           double x2 = mVPDTotEdge[i+mNVPD][ibin+1];
00416           double y1 = mVPDTotCorr[i+mNVPD][ibin];
00417           double y2 = mVPDTotCorr[i+mNVPD][ibin+1];
00418           double dcorr = y1 + (tot[i+mNVPD]-x1)*(y2-y1)/(x2-x1);
00419 
00420           mVPDLeTime[i+mNVPD] = tmp - dcorr;
00421           mTSumEast += mVPDLeTime[i+mNVPD];
00422 
00423           mVPDHitPatternEast |= 1<<i;
00424 
00425         }
00426       } else {
00427         gMessMgr->Info("","OS") << " Vpd East tube " << i+1 << " TOT out of range!" << endm;
00428         // out of range, remove this hit
00429       }
00430     }
00431   }
00432 
00433   if ( mNEast>=mVPDEastHitsCut && mNWest>=mVPDWestHitsCut ) {
00434     mVPDVtxZ = (mTSumEast/mNEast - mTSumWest/mNWest)/2.*(C_C_LIGHT/1.e9);
00435     if(fabs(mVertexZ)<200.) {
00436       mTDiff = (mTSumEast/mNEast - mTSumWest/mNWest)/2. - mVertexZ/(C_C_LIGHT/1.e9);
00437       mTStart = (mTSumEast+mTSumWest-(mNEast-mNWest)*mVertexZ/(C_C_LIGHT/1.e9))/(mNEast+mNWest);;
00438     } else {
00439       LOG_INFO << " StVpdAnalysisMaker -- vtx z " << mVertexZ << " is out of range! " << endm;
00440       mTDiff = -9999.;
00441       mTStart = -9999.;
00442     }
00443   }
00444   
00445   return;
00446 }