00001 #include "StTofMuDstReader.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 "StTriggerIdCollection.h"
00023 #include "StTriggerId.h"
00024
00025 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00026 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00027 #include "StMuDSTMaker/COMMON/StMuDst.h"
00028 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00029 #include "StMuDSTMaker/COMMON/StMuTofHit.h"
00030
00031 #include "StTofMuDstEval.h"
00032 #include "tofHitTree.h"
00033
00034
00035 ClassImp(StTofMuDstReader)
00036
00037
00038 StTofMuDstReader::StTofMuDstReader(const char *name, const char *file, StMuDstMaker* maker):StMaker(name)
00039 {
00040
00041
00042 mMuDstMaker = maker;
00043 mOutputFile = new TFile(file,"RECREATE");
00044 mOutputFile->SetCompressionLevel(1);
00045 int BufSize=(int)pow(2.,16.);
00046 int Split=1;
00047 mTree = new TTree("T","T",BufSize);
00048 mTree->SetAutoSave(10000);
00049 mTree->Branch("runId",&mTofHit.runId,"runId/I");
00050 mTree->Branch("evtId",&mTofHit.evtId,"evtId/I");
00051 mTree->Branch("vz",&mTofHit.vz,"vz/F");
00052 mTree->Branch("mult",&mTofHit.mult,"mult/I");
00053 mTree->Branch("neast",&mTofHit.neast,"neast/I");
00054 mTree->Branch("nwest",&mTofHit.nwest,"nwest/I");
00055 mTree->Branch("module",&mTofHit.module,"module/I");
00056 mTree->Branch("cell",&mTofHit.cell,"cell/I");
00057 mTree->Branch("tot",&mTofHit.tot,"tot/F");
00058 mTree->Branch("tof",&mTofHit.tof,"tof/F");
00059 mTree->Branch("L",&mTofHit.L,"L/F");
00060 mTree->Branch("beta",&mTofHit.beta,"beta/F");
00061 mTree->Branch("localY",&mTofHit.localY,"localY/F");
00062 mTree->Branch("localZ",&mTofHit.localZ,"localZ/F");
00063 mTree->Branch("deltaY",&mTofHit.deltaY,"deltaY/F");
00064 mTree->Branch("m2",&mTofHit.m2,"m2/F");
00065 mTree->Branch("trkId",&mTofHit.trkId,"trkId/I");
00066 mTree->Branch("nHitsFit",&mTofHit.nHitsFit,"nHitsFit/I");
00067 mTree->Branch("pt",&mTofHit.pt,"pt/F");
00068 mTree->Branch("eta",&mTofHit.eta,"eta/F");
00069 mTree->Branch("phi",&mTofHit.phi,"phi/F");
00070 mTree->Branch("dca",&mTofHit.dca,"dca/F");
00071 mTree->Branch("dEdx",&mTofHit.dEdx,"dEdx/F");
00072 mTree->Branch("nSigmaPi",&mTofHit.nSigmaPi,"nSigmaPi/F");
00073 mTree->Branch("nSigmaK",&mTofHit.nSigmaK,"nSigmaK/F");
00074 mTree->Branch("nSigmaP",&mTofHit.nSigmaP,"nSigmaP/F");
00075 mTree->Branch("nSigmaE",&mTofHit.nSigmaE,"nSigmaE/F");
00076 mTree->Branch("nSigmaMu",&mTofHit.nSigmaMu,"nSigmaMu/F");
00077 mTree->Branch("nSigmaTofPi",&mTofHit.nSigmaTofPi,"nSigmaTofPi/F");
00078 mTree->Branch("nSigmaTofK",&mTofHit.nSigmaTofK,"nSigmaTofK/F");
00079 mTree->Branch("nSigmaTofP",&mTofHit.nSigmaTofP,"nSigmaTofP/F");
00080 mTree->Branch("nSigmaTofE",&mTofHit.nSigmaTofE,"nSigmaTofE/F");
00081 eventNumber=0;
00082 oldEventRunId=0;
00083 }
00084
00085
00086 StTofMuDstReader::~StTofMuDstReader()
00087 {
00088 }
00089
00090
00091
00092 void StTofMuDstReader::Clear(const char*)
00093 {
00094 mTofHit.runId = -1;
00095 mTofHit.evtId = -1;
00096 mTofHit.vz = -999.;
00097 mTofHit.mult = -999;
00098 mTofHit.neast = 0;
00099 mTofHit.nwest = 0;
00100 mTofHit.module = 0;
00101 mTofHit.cell = 0;
00102 mTofHit.tot = -999.;
00103 mTofHit.tof = -999.;
00104 mTofHit.L = -999.;
00105 mTofHit.beta = -999.;
00106 mTofHit.localY = -999;
00107 mTofHit.localZ = -999.;
00108 mTofHit.deltaY = -999.;
00109 mTofHit.m2 = -999.;
00110 mTofHit.trkId = -1;
00111 mTofHit.nHitsFit = 0;
00112 mTofHit.pt = -999.;
00113 mTofHit.eta = -999.;
00114 mTofHit.phi = -999.;
00115 mTofHit.dca = -999.;
00116 mTofHit.dEdx = -999.;
00117 mTofHit.nSigmaPi = -999.;
00118 mTofHit.nSigmaK = -999.;
00119 mTofHit.nSigmaP = -999.;
00120 mTofHit.nSigmaE = -999.;
00121 mTofHit.nSigmaMu = -999.;
00122 mTofHit.nSigmaTofPi = -999.;
00123 mTofHit.nSigmaTofK = -999.;
00124 mTofHit.nSigmaTofP = -999.;
00125 mTofHit.nSigmaTofE = -999.;
00126
00127 StMaker::Clear();
00128 }
00129
00130
00131 Int_t StTofMuDstReader::Finish()
00132 {
00133 mOutputFile->Write();
00134 mOutputFile->Close();
00135 return StMaker::Finish();
00136 }
00137
00138
00139
00140 Int_t StTofMuDstReader::Init()
00141 {
00142 Clear("C");
00143 return kStOK;
00144 }
00145
00146 Int_t StTofMuDstReader::Make()
00147 {
00148 cout << "StTofMuDstReader::Make()" << endl;
00149
00150 if( mMuDstMaker != NULL ) {
00151 mMuDst = mMuDstMaker->muDst();
00152 }
00153
00154 if(!mMuDst) return kStOK;
00155
00156 StTofMuDstEval *tofEval = (StTofMuDstEval *)GetMaker("tofMuDstEval");
00157 if(!tofEval) return kStOK;
00158
00159 StMuPrimaryVertex *vtx = mMuDst->primaryVertex();
00160 if(!vtx) return kStOK;
00161
00162 StMuEvent *muEvent = mMuDst->event();
00163 if(!muEvent) return kStOK;
00164
00165 StThreeVectorF vertexPos = muEvent->primaryVertexPosition();
00166 if(fabs(vertexPos.x())<1.0e-5 && fabs(vertexPos.y())<1.0e-5 && fabs(vertexPos.z())<1.0e-5) return kStOK;
00167
00168 int runId = muEvent->runId();
00169 int evtId = muEvent->eventId();
00170 int refMult = muEvent->refMultPos() + muEvent->refMultNeg();
00171
00172 int nEast = 0;
00173 int nWest = 0;
00174 tofEval->GetPvpdNHits(nEast , nWest);
00175 cout << " nEast = " << nEast << " nWest = " << nWest << endl;
00176
00177 int nTofHits = mMuDst->numberOfTofHit();
00178 int n = mMuDst->numberOfPrimaryTracks();
00179 cout << " There are " << nTofHits << " TofHits in this event." << endl;
00180 for(int i=0;i<nTofHits;i++) {
00181 StMuTofHit *tof = (StMuTofHit*)mMuDst->tofHit(i);
00182 if(!tof) continue;
00183 int trkId = tof->associatedTrackId();
00184 int index = -1;
00185 for(int j=0;j<n;j++) {
00186 StMuTrack* t = (StMuTrack*)mMuDst->primaryTracks(j);
00187 if(!t) continue;
00188 if(trkId==t->id()) {
00189 index = j;
00190 break;
00191 }
00192 }
00193 if(index<0||index>=n) {
00194 cout << "Warning! no matched TPC track, strange!" << endl;
00195 continue;
00196 }
00197 StMuTrack *trk = (StMuTrack*)mMuDst->primaryTracks(index);
00198 if(!trk) continue;
00199
00200 Double_t local[3];
00201 tofEval->GetLocalHitPosition(tof, &local[0]);
00202 Double_t yCenter = (tof->cellIndex()-1-2.5)*3.45;
00203
00204 float m2 = pow(trk->pt()*TMath::CosH(trk->eta()),2.0)*(1./tof->beta()/tof->beta()-1.);
00205
00206 mTofHit.runId = runId;
00207 mTofHit.evtId = evtId;
00208 mTofHit.vz = vertexPos.z();
00209 mTofHit.mult = refMult;
00210 mTofHit.neast = nEast;
00211 mTofHit.nwest = nWest;
00212 mTofHit.module = tof->moduleIndex();
00213 mTofHit.cell = tof->cellIndex();
00214 mTofHit.tot = tofEval->GetUncorrectedTot(tof);
00215 mTofHit.tof = tof->timeOfFlight();
00216 mTofHit.L = tof->pathLength();
00217 mTofHit.beta = tof->beta();
00218 mTofHit.localY = local[1];
00219 mTofHit.localZ = local[2];
00220 mTofHit.deltaY = local[1]-yCenter;
00221 mTofHit.m2 = m2;
00222 mTofHit.trkId = index;
00223 mTofHit.nHitsFit = trk->nHitsFit(kTpcId)*trk->charge();
00224 mTofHit.pt = trk->pt();
00225 mTofHit.eta = trk->eta();
00226 mTofHit.phi = trk->phi();
00227 mTofHit.dca = abs(trk->dcaGlobal());
00228 mTofHit.dEdx = trk->dEdx()*1.e+6;
00229 mTofHit.nSigmaPi = trk->nSigmaPion();
00230 mTofHit.nSigmaK = trk->nSigmaKaon();
00231 mTofHit.nSigmaP = trk->nSigmaProton();
00232 mTofHit.nSigmaE = trk->nSigmaElectron();
00233 mTofHit.nSigmaMu = -999.;
00234 mTofHit.nSigmaTofPi = tof->sigmaPion()/2.;
00235 mTofHit.nSigmaTofK = tof->sigmaKaon()/2.;
00236 mTofHit.nSigmaTofP = tof->sigmaProton()/2.;
00237 mTofHit.nSigmaTofE = tof->sigmaElectron()/2.;
00238
00239 mTree->Fill();
00240
00241 }
00242
00243
00244
00245 return kStOK;
00246 }
