StRoot: StTofPool/StTofMuDstReader/StTofMuDstEval.cxx Source File

00001 #include "StTofMuDstEval.h"
00002 
00003 #include "StChain.h"
00004 #include "TF1.h"
00005 #include "TRandom.h"
00006 
00007 #include <iostream>
00008 #include <fstream>
00009 #include <stdlib.h>
00010 #include <string>
00011 #include <vector>
00012 #include <math.h>
00013 #include <cmath>
00014 
00015 #include "PhysicalConstants.h"
00016 #include "SystemOfUnits.h"
00017 #include "StThreeVector.hh"
00018 #include "StThreeVectorF.hh"
00019 #include "StThreeVectorD.hh"
00020 #include "StLorentzVectorD.hh"
00021 
00022 #include "StEvent.h"
00023 #include "StTofCollection.h"
00024 #include "StTofRawData.h"
00025 #include "StTofUtil/StSortTofRawData.h"
00026 #include "StTofUtil/StTofrGeometry.h"
00027 #include "StTriggerIdCollection.h"
00028 #include "StTriggerId.h"
00029 
00030 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00031 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00032 #include "StMuDSTMaker/COMMON/StMuDst.h"
00033 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00034 #include "StMuDSTMaker/COMMON/StMuTofHit.h"
00035 
00036 #include "StTofMuDstEval.h"
00037 
00038 ClassImp(StTofMuDstEval)
00039 
00040 //_________________________________________________
00041 StTofMuDstEval::StTofMuDstEval(const char *name, StMuDstMaker *maker):StMaker(name)
00042 {
00043   //  StTofMuDstEval Constructor
00044   // - zero all pointers defined in the header file
00045   mEvent = 0;
00046   mSortTofRawData = 0;
00047   mTofrGeom = 0;
00048   mTofCollection = 0;
00049   mMuDstMaker = maker;
00050   Clear("C");
00051 }
00052 
00053 //_________________________________________________
00054 StTofMuDstEval::~StTofMuDstEval()
00055 {
00056   Clear("C");
00057 }
00058 
00059 //_____________________________________________________________________________
00060 
00061 void StTofMuDstEval::Clear(const char*)
00062 {
00063   if(mEvent) delete mEvent;
00064   mEvent = 0;
00065   if(mSortTofRawData) delete mSortTofRawData;
00066   mSortTofRawData = 0;
00067   StMaker::Clear();
00068 }
00069 
00070 //_________________________________________________
00071 Int_t StTofMuDstEval::Finish()
00072 {
00073   if(mTofrGeom) delete mTofrGeom;
00074 
00075   return StMaker::Finish();
00076 }
00077 
00078 
00079 //_________________________________________________
00080 Int_t StTofMuDstEval::Init()
00081 {
00082   Clear("C");
00083 
00084   mTofrGeom = new StTofrGeometry("tofrGeom","tofrGeom in MuDstEval");
00085   if(!mTofrGeom->IsInitDone()) {
00086     gMessMgr->Info("TofrGemetry initialization..." ,"OS");
00087     TVolume *starHall = (TVolume *)GetDataSet("HALL");
00088     mTofrGeom->Init(starHall);
00089   }
00090 
00091   return kStOK;
00092 }
00093 //_________________________________________________
00094 Int_t StTofMuDstEval::Make()
00095 {
00096   cout << " StTofMuDstEval::Make()" << endl;
00097 
00098   if( mMuDstMaker != NULL ) {
00099     mMuDst = mMuDstMaker->muDst();
00100   }
00101 
00102   if(!mMuDst) return kStOK;
00103 
00104   mEvent = new StEvent();
00105   mTofCollection = new StTofCollection();
00106   mEvent->setTofCollection(mTofCollection);
00107   int nTofRawData = mMuDst->numberOfTofRawData();
00108   for(int i=0;i<nTofRawData;i++) {
00109     StTofRawData *aRawData;
00110     StTofRawData *aMuRawData = (StTofRawData *)mMuDst->tofRawData(i);
00111     if(aMuRawData) {
00112       unsigned short leteFlag = aMuRawData->leteFlag();
00113       unsigned short channel = aMuRawData->channel();
00114       unsigned int tdc = aMuRawData->tdc();
00115       unsigned short quality = aMuRawData->quality();
00116       aRawData = new StTofRawData(leteFlag,channel,tdc,quality);
00117     } else {
00118       aRawData = new StTofRawData(0, 0, 0, 0);
00119     }
00120     mTofCollection->addRawData(aRawData);
00121   }
00122   mSortTofRawData = new StSortTofRawData(mTofCollection);
00123 
00124   return kStOK;
00125 }
00126 //_________________________________________________
00127 void StTofMuDstEval::GetPvpdNHits(int &neast, int &nwest)
00128 {
00129   IntVec validchannel = mSortTofRawData->GetValidChannel();
00130 
00131   int used[mNPVPD]          = {0,0,0,0,0,0};
00132   int channum[mNPVPD]       ={-1,-1,-1,-1,-1,-1};
00133   float mPVPDLeTime[mNPVPD]   = {0.,0.,0.,0.,0.,0.};
00134   float mPVPDTot[mNPVPD] = {0.,0.,0.,0.,0.,0.};
00135   
00136   const Float_t mPVPDTotMin = 2.0;
00137   const Float_t mPVPDTotMax = 8.0;
00138 
00139   for(unsigned int ich=0;ich<validchannel.size();ich++){
00140     int chan = validchannel[ich];
00141     if(chan<mNTOFr5) continue;       // need only pvpd 
00142     int ichan = chan - mNTOFr5;
00143     if(ichan<0||ichan>=mNPVPD) continue;
00144     if(used[ichan]>0) continue;                   // skip multi hits
00145     used[ichan]++;
00146     // leading edge
00147     int tmptdc = (mSortTofRawData->GetLeadingTdc(chan))[0];    
00148     float pvpdletdc = tmptdc;
00149     mPVPDLeTime[ichan] = pvpdletdc * VHRBIN2PS/1000.;
00150     // Trailing edge
00151     tmptdc = (mSortTofRawData->GetTrailingTdc(chan))[0];    
00152     float pvpdtetdc = tmptdc;
00153     float tetime = pvpdtetdc * HRBIN2PS/1000.;
00154     mPVPDTot[ichan] = tetime - mPVPDLeTime[ichan];
00155     channum[ichan]=ichan;
00156   }
00157   
00158   // sort out the piled-up hit - remove those hits not from this event
00159   int nPVPDFired = 0;
00160   for(int i=0;i<mNPVPD;i++) {
00161     if(channum[i]!=i) {
00162       mPVPDLeTime[i] = 0.;
00163       mPVPDTot[i] = 0.;
00164     } else {
00165       nPVPDFired++;
00166     }
00167   }
00168   for(int i=0;i<mNPVPD;i++) {
00169     if(channum[i]!=i) continue;
00170     int n0 = 0;
00171     for(int j=0;j<mNPVPD;j++) {
00172       if(channum[j]!=j) continue;
00173       float dt = mPVPDLeTime[j] - mPVPDLeTime[i];
00174       if(fabs(dt)<200.) n0++;
00175     }
00176 
00177     if(n0>=nPVPDFired/2) {  // OK
00178     } else { // not from this event
00179       mPVPDLeTime[i] = 0.;
00180       mPVPDTot[i] = 0.;
00181     }
00182   }
00183 
00184   
00185   for(int i=0;i<mNPVPD/2;i++) {
00186     if(mPVPDLeTime[i]>0.&&mPVPDTot[i]>mPVPDTotMin&&mPVPDTot[i]<mPVPDTotMax) neast++;
00187     int j=i+mNPVPD/2;
00188     if(mPVPDLeTime[j]>0.&&mPVPDTot[j]>mPVPDTotMin&&mPVPDTot[j]<mPVPDTotMax) nwest++;
00189   }
00190 
00191   return;
00192 }
00193 //_________________________________________________
00194 Float_t StTofMuDstEval::GetUncorrectedTot(StMuTofHit *tofHit)
00195 {
00196   Float_t tot = -999.;
00197 
00198   int daqIndex = tofHit->daqIndex();
00199   IntVec validchannel = mSortTofRawData->GetValidChannel();
00200   for(unsigned int ich=0;ich<validchannel.size();ich++){
00201     int chan = validchannel[ich];
00202     if(chan==daqIndex) {
00203       int tmptdc = (mSortTofRawData->GetLeadingTdc(chan))[0];    
00204       float letime = tmptdc * VHRBIN2PS/1000.;
00205       tmptdc = (mSortTofRawData->GetTrailingTdc(chan))[0];    
00206       float tetime = tmptdc * HRBIN2PS/1000.;
00207       tot = tetime - letime;
00208       break;
00209     }
00210   }
00211 
00212   return tot;
00213 }
00214 //_________________________________________________
00215 void StTofMuDstEval::GetLocalHitPosition(StMuTofHit *tofHit, Double_t* local)
00216 {
00217   if(!tofHit) return;
00218   int itray = tofHit->trayIndex();
00219   int imodule = tofHit->moduleIndex();
00220   StTofrGeomSensor* sensor = mTofrGeom->GetGeomSensor(imodule,itray);
00221   if(!sensor) {
00222     gMessMgr->Warning("","OS") << " No sensitive module in the projection??? -- Something weird!!! " << endm;
00223     return;
00224   }
00225 
00226   Double_t global[3] = {0.,0.,0.};
00227   global[0] = tofHit->projectedPoint().x();
00228   global[1] = tofHit->projectedPoint().y();
00229   global[2] = tofHit->projectedPoint().z();
00230 
00231   sensor->Master2Local(&global[0],local);
00232   delete sensor;
00233   return;
00234 }