StRoot: StVpdCalibMaker/StVpdCalibMaker.cxx Source File

00001 /*******************************************************************
00002  *
00003  * $Id: StVpdCalibMaker.cxx,v 1.11 2011/02/23 20:00:52 geurts Exp $
00004  *
00005  * Author: Xin Dong
00006  *****************************************************************
00007  *
00008  * Description: - VPD Calibration Maker to do the calibration for upVPD 
00009  *              - store into StBTofHeader
00010  *
00011  *****************************************************************
00012  *
00013  * $Log: StVpdCalibMaker.cxx,v $
00014  * Revision 1.11  2011/02/23 20:00:52  geurts
00015  * Change MaxBin for ToT arrays from 60 to 128 (in agreement with the IDL definition of vpdTotCorr)
00016  * Move the log message that informs the user about the start-timing mode outside the tube loop ... no need to see the same message 38 times.
00017  *
00018  * Revision 1.10  2011/01/07 21:28:51  geurts
00019  * Allow user to "force" VPD-start or startless mode, regardless of dbase setting
00020  *
00021  * Revision 1.9  2010/12/18 01:10:27  geurts
00022  * bugfix: apply VPD outlier correction in StEvent mode, too.
00023  *
00024  * Revision 1.8  2010/12/10 21:38:06  geurts
00025  * RT#1996: default initialization of database table in case dataset is not found
00026  *
00027  * Revision 1.7  2010/10/31 05:43:57  geurts
00028  * apply outlier truncation to all energies, keep at 20% (previously only enabled for low beam energies)
00029  *
00030  * Revision 1.6  2010/05/25 22:09:48  geurts
00031  * improved database handling and reduced log output
00032  *
00033  * Revision 1.5  2010/05/12 22:46:51  geurts
00034  * Startless BTOF self-calibration method (Xin)
00035  *
00036  * Revision 1.4  2010/05/06 22:37:41  geurts
00037  * Remove slower hits (outliers) in VPD timing calculations (Xin Dong)
00038  *
00039  * Revision 1.3  2010/01/29 19:51:08  geurts
00040  * Fixed a bug in vzVpdFinder outlier bookkeeping (Xin Dong)
00041  *
00042  * Revision 1.2  2009/12/04 22:22:17  geurts
00043  * two updates (Xin):
00044  * - use the new calibration table vpdTotCorr
00045  * - update on the vzVpdFinder(), a cut on timing diff is used to remove outliers.
00046  *
00047  * Revision 1.1  2009/11/09 20:55:54  geurts
00048  * first release
00049  *
00050  *
00051  *******************************************************************/
00052 #include <iostream>
00053 #include "TFile.h"
00054 
00055 #include "StEvent.h"
00056 #include "StBTofCollection.h"
00057 #include "StBTofHit.h"
00058 #include "StBTofHeader.h"
00059 #include "StEventTypes.h"
00060 #include "Stypes.h"
00061 #include "StMessMgr.h"
00062 #include "StThreeVectorD.hh"
00063 #include "StEventUtilities/StuRefMult.hh"
00064 #include "PhysicalConstants.h"
00065 #include "phys_constants.h"
00066 #include "tables/St_vpdTotCorr_Table.h"
00067 
00068 #include "StBTofUtil/StBTofHitCollection.h"
00069 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00070 #include "StMuDSTMaker/COMMON/StMuDst.h"
00071 #include "StMuDSTMaker/COMMON/StMuBTofHit.h"
00072 
00073 #include "StVpdCalibMaker.h"
00074 
00075 #include "StMessMgr.h"
00076 #include "StMemoryInfo.hh"
00077 #include "StTimer.hh"
00078 
00079 #ifdef __ROOT__
00080 ClassImp(StVpdCalibMaker)
00081 #endif
00082 
00084   const Double_t StVpdCalibMaker::VHRBIN2PS =  24.4140625; // 1000*25/1024 (ps/chn)
00086   const Double_t StVpdCalibMaker::HRBIN2PS = 97.65625; // 97.65625= 1000*100/1024  (ps/chn)
00088   const Double_t StVpdCalibMaker::TMAX = 51200.;
00090   const Double_t StVpdCalibMaker::VZDIFFCUT=6.;
00092   const Double_t StVpdCalibMaker::TDIFFCUT=0.8;
00094   const Double_t StVpdCalibMaker::FracTruncated=0.2;
00095 
00096 //_____________________________________________________________________________
00097 StVpdCalibMaker::StVpdCalibMaker(const Char_t *name) : StMaker(name)
00098 {
00102   setVPDHitsCut(1,1);
00103 
00104   //mEvent   = 0;
00105   mTruncation = kFALSE;  
00106   mMuDst   = 0;
00107   mMuDstIn = kFALSE;   
00108   mBTofColl       = 0;
00109   mValidCalibPar  = kFALSE;
00110 
00111   setCreateHistoFlag(kFALSE);
00112   setHistoFileName("vpdana.root");
00113 
00114   // default initialization from database
00115   mInitFromFile = kFALSE;
00116   // assign default locations and names to the calibration files 
00117   setCalibFilePvpd("/star/institutions/rice/calib/default/pvpdCali_4DB.dat");
00118   // use vpd as start by default;
00119   mUseVpdStart = kTRUE;
00120   mForceTofStart = kFALSE; // flag indicates user-override for TOF Start time calculation
00121 }
00122 
00123 //_____________________________________________________________________________
00124 StVpdCalibMaker::~StVpdCalibMaker()
00125 { /* noop */ }
00126 
00127 //_____________________________________________________________________________
00128 void StVpdCalibMaker::resetPars()
00129 {
00130   memset(mVPDTotEdge, 0, sizeof(mVPDTotEdge));
00131   memset(mVPDTotCorr, 0, sizeof(mVPDTotCorr));
00132 }
00133 
00134 //_____________________________________________________________________________
00135 void StVpdCalibMaker::resetVpd()
00136 {
00137   memset(mVPDLeTime, 0, sizeof(mVPDLeTime));
00138   memset(mVPDTot, 0, sizeof(mVPDTot));
00139   memset(mFlag, 0, sizeof(mFlag));
00140 
00141   for(int i=0;i<MaxVpdVz;i++) {
00142     mVPDVtxZ[i] = -9999.;
00143   }
00144   mNVzVpd = 0;
00145   mVPDHitPatternEast = 0;
00146   mVPDHitPatternWest = 0;
00147 }
00148 
00149 //____________________________________________________________________________
00150 Int_t StVpdCalibMaker::Init()
00151 {
00152   resetPars();
00153   resetVpd();
00154 
00155   // m_Mode can be set by SetMode() method
00156   if(m_Mode) {
00157 //    setHistoFileName("vpdana.root");
00158   } else {
00159     setHistoFileName("");    
00160   }
00161 
00162   if (mHisto){
00163     bookHistograms();
00164     LOG_INFO << "Histograms are booked" << endm;
00165     if (mHistoFileName!="") {
00166       LOG_INFO << "Histograms will be stored in " << mHistoFileName.c_str() << endm;
00167     }
00168   }
00169 
00170   //return kStOK;
00171   return StMaker::Init();
00172 }
00173 
00174 //____________________________________________________________________________
00175 Int_t StVpdCalibMaker::InitRun(Int_t runnumber)
00176 {
00177   // tof run configurations
00178 
00179   Int_t val = kStOK;
00180   val = initParameters(runnumber);
00181   if(val==kStOK) {
00182     mValidCalibPar = kTRUE;
00183     //return kStOK;
00184     return StMaker::InitRun(runnumber);
00185   } else {
00186     mValidCalibPar = kFALSE;
00187     LOG_WARN << " !!! No valid calibration parameters for VPD! " << endm;
00188     return kStWarn;
00189   }
00190 
00191 
00192 }
00193 
00194 //_____________________________________________________________________________
00195 Int_t StVpdCalibMaker::initParameters(Int_t runnumber)
00196 {
00199   if (mInitFromFile){
00200     LOG_INFO << "Initializing VPD calibration parameters from file"
00201                  << "(" << mCalibFilePvpd << ")" << endm;  
00202     ifstream inData;
00203     inData.open(mCalibFilePvpd.c_str());
00204     int nchl, nbin;
00205     for(int i=0;i<NVPD*2;i++) {
00206       inData>>nchl;
00207       inData>>nbin;
00208       if (nbin>NBinMax) {
00209                 LOG_ERROR << "nummer of bins (" << nbin << ") out of range ("
00210                           << NBinMax << ") for vpd channel " << i << endm;
00211                 return kStErr;
00212       } 
00213       for(int j=0;j<=nbin;j++) inData>>mVPDTotEdge[i][j];
00214       for(int j=0;j<=nbin;j++) inData>>mVPDTotCorr[i][j];
00215     }
00216     inData.close();
00217   }
00218   else {
00220     LOG_INFO << "Initializing VPD calibration parameters from database" << endm;
00221 
00222 // read in and check the size
00223 //     TDataSet *dbDataSet = GetDataBase("Calibrations/tof");
00224 //     if (!dbDataSet){
00225 //       LOG_ERROR  << "unable to get VPD/TOF run parameters" << endm;
00226 //       return kStErr;
00227 //     }
00228   
00229     // read vpdTotCorr table
00230     TDataSet *dbDataSet = GetDataBase("Calibrations/tof/vpdTotCorr");
00231     if (dbDataSet){
00232       St_vpdTotCorr* vpdTotCorr = static_cast<St_vpdTotCorr*>(dbDataSet->Find("vpdTotCorr"));
00233       if(!vpdTotCorr) {
00234         LOG_ERROR << "unable to get vpdTotCorr table parameters" << endm;
00235         //    assert(vpdTotCorr);
00236         return kStErr;
00237       }
00238       vpdTotCorr_st* totCorr = static_cast<vpdTotCorr_st*>(vpdTotCorr->GetArray());
00239       Int_t numRows = vpdTotCorr->GetNRows();
00240 
00241       if(numRows!=NVPD*2) {
00242         LOG_WARN  << " Mis-matched number of rows in vpdTotCorr table: " << numRows 
00243                   << " (exp:" << NVPD*2 << ")" << endm;
00244       }
00245 
00246       LOG_DEBUG << " Number of rows read in: " << numRows << " for Vpd ToT correction" << endm;
00247 
00248       for (Int_t i=0;i<numRows;i++) {
00249         if (!mForceTofStart){
00250           if(i==0) {  // identify once only, for the first tube
00251             short flag = totCorr[i].corralgo;
00252             if(flag==0) {
00253               mUseVpdStart=kTRUE;
00254               LOG_INFO << "Selected VPD for TOF start-timing (corralgo=1)" << endm;
00255             }
00256             else if(flag==1) {
00257               mUseVpdStart=kFALSE;
00258               LOG_INFO << "VPD NOT used for TOF start-timing (corralgo=0)" << endm;
00259             }
00260             else {
00261               LOG_WARN << "Unknown calibration option " << flag << endm;
00262             }
00263           }
00264           else { // verify that all other entries agree
00265             if (totCorr[0].corralgo==0 && (totCorr[i].corralgo!=0))
00266               {LOG_WARN << "corralgo dbase inconsistency: " << totCorr[i].corralgo << endm;}
00267             if (totCorr[0].corralgo==1 && (totCorr[i].corralgo!=1))
00268               {LOG_WARN << "corralgo dbase inconsistency: " << totCorr[i].corralgo << endm;}
00269           }
00270         }
00271 
00272         short tubeId = totCorr[i].tubeId;
00273         // check index range
00274         if (tubeId>2*NVPD) {
00275           LOG_ERROR << "tubeId (" << tubeId << ") out of range ("
00276                     << 2*NVPD << ")" << endm;
00277           return kStErr;
00278         }
00279 
00280         for(Int_t j=0;j<NBinMax;j++) {
00281           mVPDTotEdge[tubeId-1][j] = totCorr[i].tot[j];
00282           mVPDTotCorr[tubeId-1][j] = totCorr[i].corr[j];
00283           LOG_DEBUG << " east/west: " << (tubeId-1)/NVPD << " tubeId: " << tubeId << endm;
00284         } // end j 0->NBinMax
00285       } // end i 0->numRows
00286 
00287       // Let the user know what is used for calculating TOF Start (and how we got there)
00288       if (mForceTofStart) { LOG_INFO << "Detected user override in Start Timing selection ::setUseVpdStart()." << endm;}
00289       if (mUseVpdStart) {
00290         LOG_INFO << "Selected VPD for TOF start-timing" << endm;
00291       }
00292       else {
00293         LOG_INFO << "VPD NOT used for TOF start-timing" << endm;  
00294       }
00295     }
00296     else {
00297     // Note: this construction addresses RT#1996 and allows backward compatibility with older database
00298     //       timestamps at which these database structures did not exist. All values will be zero (hence the ERROR)
00299     //       and the VPD will be disabled for TOF start timing, i.e. the TOF could still operate in start-less mode 
00300       LOG_ERROR << "unable to get vpdTotCorr dataset ... reset all to zero values (NOT GOOD!) and disable use for TOF-start" << endm; 
00301       resetPars();
00302       mUseVpdStart=kFALSE;
00303       LOG_INFO << "VPD NOT used for TOF start-timing" << endm;  
00304       return kStErr;
00305     }
00306   }
00307 
00308   return kStOK;
00309 }
00310 
00311 
00312 //_____________________________________________________________________________
00313 Int_t StVpdCalibMaker::Finish()
00314 {
00315   if (mHistoFileName!="") writeHistograms();
00316   return StMaker::Finish();//kStOK;
00317 }
00318 
00319 
00320 //_____________________________________________________________________________
00321 Int_t StVpdCalibMaker::Make()
00322 {
00323   LOG_DEBUG << " StVpdCalibMaker::Make: starting ..." << endm;
00324   if(mHisto&&mhEventCounter) mhEventCounter->Fill(0);
00325   if(!mValidCalibPar) {
00326     LOG_WARN << " No valid calibration parameters. Skip ... " << endm;
00327     return kStWarn;
00328   }
00329 
00330   if(mHisto&&mhEventCounter) mhEventCounter->Fill(1);
00331   resetVpd();
00332   bool loadOK = loadVpdData();
00333 
00334   if(!loadOK) return kStWarn;
00335   if(mHisto&&mhEventCounter) mhEventCounter->Fill(2);
00336 
00337   tsum(mVPDTot, mVPDLeTime);
00338   vzVpdFinder();
00339 
00340   if(mHisto) fillHistograms();
00341   if(mHisto&&mhEventCounter) mhEventCounter->Fill(3);
00342   bool writeOK = writeVpdData();
00343   if (!writeOK) return kStWarn;
00344   if(writeOK&&mHisto) mhEventCounter->Fill(4);
00345 
00346   return StMaker::Make(); //kStOK;
00347 }
00348 
00349 
00350 //_____________________________________________________________________________
00351 Bool_t StVpdCalibMaker::writeVpdData() const
00352 {
00353   StBTofHeader *tofHeader = 0;
00354   if(mMuDstIn) {
00355     if(!mMuDst) {
00356       LOG_WARN << " No MuDst to write ... " << endm;
00357       return kFALSE;
00358     }
00359     tofHeader = (StBTofHeader *) mMuDst->btofHeader();
00360   } else {
00361     if(!mBTofColl) {
00362       LOG_WARN << " No StEvent/btofCollection to write ... " << endm;
00363       return kFALSE;
00364     }
00365     tofHeader = (StBTofHeader *) mBTofColl->tofHeader();
00366   }
00367 
00368   if(!tofHeader) return kFALSE;
00370   for(int i=0;i<NVPD;i++) {
00371     tofHeader->setVpdTime(west, i+1, mVPDLeTime[i]);
00372     tofHeader->setVpdTime(east, i+1, mVPDLeTime[i+NVPD]);
00373   }
00374   tofHeader->setVpdHitPattern(east, mVPDHitPatternEast);
00375   tofHeader->setVpdHitPattern(west, mVPDHitPatternWest);
00376   for(int i=0;i<mNVzVpd;i++) {
00377     tofHeader->setVpdVz(mVPDVtxZ[i],i);
00378   }
00379 
00380 
00381     LOG_INFO << "BTofHeader: NWest = " << tofHeader->numberOfVpdHits(west) 
00382              << " NEast = " << tofHeader->numberOfVpdHits(east) << endm;
00383     if(tofHeader->numberOfVpdHits(west)!=mNWest ||
00384        tofHeader->numberOfVpdHits(east)!=mNEast) {
00385       LOG_WARN << "BTofHeader inconsistency: Local nWest = " << mNWest << " nEast = " << mNEast << endm;
00386     }
00387     LOG_INFO <<"summary:  VPD-VtxZ = " << mVPDVtxZ[0] 
00388              << "; TSum West = " << mTSumWest << "  East = " << mTSumEast << endm;
00389 
00390   return kTRUE;
00391 }
00392 
00393 
00394 //____________________________________________________________________________
00395 Bool_t StVpdCalibMaker::loadVpdData()
00396 {
00397   if(mMuDstIn) {
00398     StMuDstMaker *mMuDstMaker = (StMuDstMaker *)GetMaker("MuDst");
00399     if(!mMuDstMaker) {
00400       LOG_WARN << "No MuDstMaker ... bye-bye" << endm;
00401       return kFALSE;
00402     }
00403     mMuDst = mMuDstMaker->muDst();
00404     if(!mMuDst) {
00405       LOG_WARN << "No MuDst ... bye-bye" << endm;
00406       return kFALSE;
00407     }
00408 
00409     // Apply truncation to all energies and collision systems
00410 //    if(mMuDst->event() && mMuDst->event()->runInfo().centerOfMassEnergy()<40.) {
00411       mTruncation = kTRUE;
00412 //    } else {
00413 //      mTruncation = kFALSE;
00414 //    }
00415 
00416 
00417     Int_t nhits = mMuDst->numberOfBTofHit();
00418 //    Int_t nhits = mMuDst->btofArray(muBTofHit)->GetEntries();
00419     for(int i=0;i<nhits;i++) {
00420       StMuBTofHit *aHit = (StMuBTofHit *)mMuDst->btofHit(i);
00421       if(!aHit) continue;
00422       int trayId = aHit->tray();
00423        if(trayId==WestVpdTrayId || trayId==EastVpdTrayId) {   // VPD this time
00424         int tubeId = aHit->cell();
00425         int indx = (trayId-NTray-1)*NVPD+(tubeId-1);
00426         if (indx>=2*NVPD){
00427           LOG_ERROR << "vpd index (" << indx << ") out of range ("
00428                     << 2*NVPD << ") for trayId-tubeId " << trayId 
00429                     << "-" << tubeId << endm;
00430           return kStErr;
00431         }
00432         mVPDLeTime[indx] = aHit->leadingEdgeTime();
00433         mVPDTot[indx] = aHit->tot();
00434       }
00435     }
00436 
00437   } else {
00438     StEvent *thisEvent = (StEvent *) GetInputDS("StEvent");
00439   
00440     // event selection  // no primary vertex required
00441     if( !thisEvent )  {
00442       LOG_WARN << "No StEvent present" << endl;
00443       return kFALSE;
00444      }
00445     if (!thisEvent->btofCollection() ) {
00446        LOG_WARN << "No BTOFCollection present" << endm;
00447        return kFALSE;
00448     }
00449     if (!thisEvent->btofCollection()->hitsPresent() ) {
00450       LOG_WARN << "No hits present" << endm;
00451       return kFALSE;
00452     }
00453 
00454 //    if(thisEvent->runInfo() && thisEvent->runInfo()->centerOfMassEnergy()<40.) { // 40 GeV or less
00455       mTruncation = kTRUE;
00456 //    } else {
00457 //      mTruncation = kFALSE;
00458 //    }
00459 
00460     mBTofColl = thisEvent->btofCollection();
00461     StSPtrVecBTofHit &tofHits = mBTofColl->tofHits();
00462     Int_t nhits = tofHits.size();
00463     LOG_INFO << "Total number of TOF cells + VPD tubes : " << nhits << endm;
00464   
00465     for(int i=0;i<nhits;i++) {
00466       StBTofHit *aHit = dynamic_cast<StBTofHit*>(tofHits[i]);
00467       if(!aHit) continue;
00468       int trayId = aHit->tray();
00469       if(trayId==WestVpdTrayId || trayId==EastVpdTrayId) {   // VPD this time
00470         int tubeId = aHit->cell();
00471         int indx = (trayId-NTray-1)*NVPD+(tubeId-1);
00472         if (indx>=2*NVPD){
00473           LOG_ERROR << "vpd index (" << indx << ") out of range ("
00474                     << 2*NVPD << ") for trayId-tubeId " << trayId 
00475                     << "-" << tubeId << endm;
00476           return kStErr;
00477         }
00478         mVPDLeTime[indx] = aHit->leadingEdgeTime();
00479         mVPDTot[indx] = aHit->tot();
00480       }
00481     } // end loop hits
00482   }
00483 
00484   return kTRUE;
00485 }
00486 
00487 
00488 //_____________________________________________________________________________
00489 void StVpdCalibMaker::tsum(const Double_t *tot, const Double_t *time)
00490 {
00492   mTSumEast = 0.;
00493   mTSumWest = 0.;
00494   mNEast = 0;
00495   mNWest = 0;
00496   mVPDHitPatternWest = 0;
00497   mVPDHitPatternEast = 0;
00498 
00499   bool vpdEast=false;
00500   for(int i=0;i<2*NVPD;i++) {
00501     if (i>=NVPD) vpdEast=true;
00502 
00503     if( time[i]>0. && tot[i]>0. ) {
00504 
00505       int ibin = -1;
00506       for(int j=0;j<NBinMax-1;j++) {
00507         if(tot[i]>=mVPDTotEdge[i][j] && tot[i]<mVPDTotEdge[i][j+1]) {
00508           ibin = j;
00509           break;
00510         }
00511       }
00512       if(ibin>=0&&ibin<NBinMax) {
00513         Double_t x1 = mVPDTotEdge[i][ibin];
00514         Double_t x2 = mVPDTotEdge[i][ibin+1];
00515         Double_t y1 = mVPDTotCorr[i][ibin];
00516         Double_t y2 = mVPDTotCorr[i][ibin+1];
00517         Double_t dcorr = y1 + (tot[i]-x1)*(y2-y1)/(x2-x1);
00518 
00519         if (vpdEast){
00520           mNEast++;
00521           mVPDLeTime[i] = time[i] - dcorr;
00522           mTSumEast += mVPDLeTime[i];
00523           mVPDHitPatternEast |= 1<<(i-NVPD);
00524           mFlag[i] = 1;
00525         } else {
00526           mNWest++;
00527           mVPDLeTime[i] = time[i] - dcorr;
00528           mTSumWest += mVPDLeTime[i];
00529           mVPDHitPatternWest |= 1<<i;
00530           mFlag[i] = 1;
00531         }
00532 
00533       } else {
00534         if (vpdEast){
00535           LOG_WARN << " Vpd East tube " << i+1-NVPD << " TOT ("<< ibin
00536                    << ") out of range (0-"<<NBinMax<<") !" << endm;
00537         } else{
00538           LOG_WARN << " Vpd West tube " << i+1 << " TOT ("<< ibin
00539                    << ") out of range (0-"<<NBinMax<<") !" << endm;
00540         }
00541         mVPDLeTime[i] = 0.; // out of range, remove this hit
00542       }
00543     }
00544   }
00545 
00546 }
00547 
00548 
00549 //_____________________________________________________________________________
00551 void StVpdCalibMaker::vzVpdFinder()
00552 {
00553   for(int i=0;i<2*NVPD;i++){
00554     // check
00555     if(mVPDLeTime[i]<1.e-4 || !mFlag[i]) continue;
00556     double vpdtime;
00557     if(i<NVPD&&mNWest>1) {  // west VPD
00558       vpdtime = (mVPDLeTime[i]*mNWest-mTSumWest)/(mNWest-1);
00559       if(fabs(vpdtime)>TDIFFCUT) {
00560         mTSumWest -= mVPDLeTime[i];
00561         mVPDLeTime[i] = 0.;
00562         mNWest--;
00563         mVPDHitPatternWest &= ( 0x7ffff - (1<<i) );
00564         mFlag[i] = 0;
00565       }
00566     }
00567     if(i>=NVPD&&mNEast>1) {  // east VPD
00568       vpdtime = (mVPDLeTime[i]*mNEast-mTSumEast)/(mNEast-1);
00569       if(fabs(vpdtime)>TDIFFCUT) {
00570         mTSumEast -= mVPDLeTime[i];
00571         mVPDLeTime[i] = 0.;
00572         mNEast--; 
00573         mVPDHitPatternEast &= ( 0x7ffff - (1<<(i-NVPD)) );
00574         mFlag[i] = 0;
00575       }
00576     }
00577   }
00578 
00579   // remove slower hit in low energy runs.
00580   if(mTruncation) {
00581     Int_t hitIndex[2*NVPD];
00582     Int_t nTube = NVPD;
00583     TMath::Sort(nTube, &mVPDLeTime[0], &hitIndex[0]);
00584     int nRejectedWest = (int)(FracTruncated*mNWest+0.5);
00585     LOG_INFO << " NWest before = " << mNWest << " rejected = " << nRejectedWest << endm;
00586     for(int i=0;i<nRejectedWest;i++) {
00587       int index = hitIndex[i];
00588       mTSumWest -= mVPDLeTime[index];
00589       mVPDLeTime[index] = 0.;
00590       mNWest--;
00591       mVPDHitPatternWest &= ( 0x7ffff - (1<<index) );
00592       mFlag[index] = 0;
00593     }
00594 
00595     TMath::Sort(nTube, &mVPDLeTime[NVPD], &hitIndex[NVPD]);
00596     int nRejectedEast = (int)(FracTruncated*mNEast+0.5);
00597     LOG_INFO << " NEast before = " << mNEast << " rejected = " << nRejectedEast << endm;
00598     for(int i=0;i<nRejectedEast;i++) {
00599       int index = hitIndex[i+NVPD] + NVPD;
00600       mTSumEast -= mVPDLeTime[index];
00601       mVPDLeTime[index] = 0.;
00602       mNEast--;
00603       mVPDHitPatternEast &= ( 0x7ffff - (1<<(index-NVPD)) );
00604       mFlag[index] = 0;
00605     }
00606   }
00607 
00608   // calculate the vertex z from vpd
00609   if ( mNEast>=mVPDEastHitsCut && mNWest>=mVPDWestHitsCut ) {
00610     mVPDVtxZ[0] = (mTSumEast/mNEast - mTSumWest/mNWest)/2.*(C_C_LIGHT/1.e9);
00611     mNVzVpd++;
00612   }
00613 }
00614 
00615 
00616 //_____________________________________________________________________________
00617 void StVpdCalibMaker::bookHistograms()
00618 {
00619   mhEventCounter = new TH1D("eventCounter","eventCounter",20,0,20);
00620   mhNVpdHits = new TH2D("vpdHits"," west vs east ",20,0.,20.,20,0.,20.);
00621   for(int i=0;i<2*NVPD;i++) {
00622     //char buf[100];
00623     TString buf;
00624     if(i<NVPD) {
00625       //sprintf(buf,"res_W_%d",i+1);
00626       buf.Form("res_W_%d",i+1);
00627     } else {
00628       //sprintf(buf,"res_E_%d",i-mNVPD+1);
00629       buf.Form("res_E_%d",i-NVPD+1);
00630     }
00631     mhVpd[i] = new TH1D(buf, buf, 3000, -3., 3.);
00632   }
00633   mhVpdAll = new TH1D("res_All","res_All", 3000, -3., 3.);
00634 }
00635 
00636 //_____________________________________________________________________________
00637 void StVpdCalibMaker::fillHistograms()
00638 {
00639   if(!mHisto) return;
00640   if (mhNVpdHits) mhNVpdHits->Fill(mNEast, mNWest);
00641   if(mNWest>=2) {
00642     for(int i=0;i<NVPD;i++) {
00643       if(mVPDLeTime[i]>0.) {
00644         Double_t tdiff = (mNWest*mVPDLeTime[i] - mTSumWest)/(mNWest - 1);
00645         if (mhVpd[i]) mhVpd[i]->Fill(tdiff);
00646         if (mhVpdAll) mhVpdAll->Fill(tdiff);
00647       }
00648     }
00649   }
00650   if(mNEast>=2) {
00651     for(int j=0;j<NVPD;j++) {
00652       int i = j + NVPD;
00653       if(mVPDLeTime[i]>0.) {
00654         Double_t tdiff = (mNWest*mVPDLeTime[i] - mTSumEast)/(mNEast - 1);
00655         if (mhVpd[i]) mhVpd[i]->Fill(tdiff);
00656         if (mhVpdAll) mhVpdAll->Fill(tdiff);
00657       }
00658     }
00659   }
00660 }
00661 
00662 
00663 //_____________________________________________________________________________
00664 void StVpdCalibMaker::writeHistograms() const
00665 {
00666   // Output file
00667   TFile *theHistoFile =  new TFile(mHistoFileName.c_str(), "RECREATE");
00668   LOG_INFO << "StVpdCalibMaker::writeHistograms()"
00669        << " histogram file " <<  mHistoFileName << endm;
00670 
00671   if (theHistoFile&&theHistoFile->IsOpen()){
00672     theHistoFile->cd();
00673 
00674     if(mHisto) {
00675       mhEventCounter->Write();
00676       mhNVpdHits->Write();
00677       mhVpdAll->Write();
00678       for(int i=0;i<2*NVPD;i++) {
00679         mhVpd[i]->Write();
00680       }
00681     }
00682 
00683     theHistoFile->Close();
00684     delete theHistoFile;
00685   }
00686   else 
00687     LOG_ERROR << "unable to open histogram file" << endm;
00688     
00689   return;
00690 }