StRoot: StJetMaker/emulator/StjeDefaultJetTreeWriter.cxx Source File

00001 // $Id: StjeDefaultJetTreeWriter.cxx,v 1.12 2010/11/29 02:40:26 pibero Exp $
00002 // Copyright (C) 2008 Tai Sakuma <sakuma@bnl.gov>
00003 #include "StjeDefaultJetTreeWriter.h"
00004 
00005 #include <StMuDSTMaker/COMMON/StMuDst.h>
00006 #include <StMuDSTMaker/COMMON/StMuEvent.h>
00007 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
00008 
00009 #include "StSpinPool/StJets/StJets.h"
00010 #include "StSpinPool/StJets/StJet.h"
00011 #include "StMuTrackFourVec.h"
00012 
00013 #include "StMuTrackEmu.h"
00014 #include "StMuTowerEmu.h"
00015 #include "StMcTrackEmu.h"
00016 
00017 #include <StJetFinder/StProtoJet.h>
00018 
00019 
00020 #include <StFourPMaker.h>
00021 #include <StBET4pMaker.h>
00022 
00023 #include <TTree.h>
00024 
00025 using namespace std;
00026 
00027 StjeDefaultJetTreeWriter::StjeDefaultJetTreeWriter(StMuDstMaker& uDstMaker, std::string outFileName)
00028   : _uDstMaker(uDstMaker)
00029   , _OutFileName(outFileName)
00030   , _jetTree(0)
00031   , _outFile(0)
00032 {
00033 
00034 }
00035 
00036 StjeDefaultJetTreeWriter::~StjeDefaultJetTreeWriter()
00037 {
00038 
00039 }
00040 
00041 void StjeDefaultJetTreeWriter::addJetFinder(StFourPMaker* fourPMaker, const vector<const AbstractFourVec*>* particleList, list<StProtoJet>* protoJetList, const char* name, StJets* stjets)
00042 {
00043   AnalyzerCtl anaCtl;
00044   anaCtl._branchName = name;
00045   anaCtl._fourPMaker = fourPMaker;
00046   anaCtl._protoJetList = protoJetList;
00047   anaCtl._jets = stjets;
00048 
00049   _analyzerCtlList.push_back(anaCtl);
00050 }
00051 
00052 void StjeDefaultJetTreeWriter::Init()
00053 {
00054   _outFile = new TFile(_OutFileName.c_str(), "recreate");
00055   _jetTree  = new TTree("jet", "jetTree");
00056 
00057   for(vector<AnalyzerCtl>::iterator it = _analyzerCtlList.begin(); it != _analyzerCtlList.end(); ++it) {
00058     _jetTree->Branch((*it)._branchName.c_str(), "StJets", &((*it)._jets));
00059   }
00060 
00061   _jetTree->BranchRef();
00062 }
00063 
00064 void StjeDefaultJetTreeWriter::Finish()
00065 {
00066   _outFile->Write();
00067   _outFile->Close();
00068   delete _outFile;
00069   _outFile = 0;
00070 }
00071 
00072 void StjeDefaultJetTreeWriter::fillJetTreeHeader(int iAnalyzer)
00073 {
00074   StFourPMaker* fourPMaker = _analyzerCtlList[iAnalyzer]._fourPMaker;
00075   StJets& jets = *_analyzerCtlList[iAnalyzer]._jets;
00076 
00077   jets.Clear();
00078   jets.setBemcCorrupt(fourPMaker->bemcCorrupt());
00079 
00080   StMuEvent* event = _uDstMaker.muDst()->event();
00081   jets.seteventId(event->eventId());
00082   jets.seteventNumber(event->eventNumber());
00083   jets.setrunId(event->runId());
00084   jets.setrunNumber(event->runNumber());
00085 
00086   StBET4pMaker* bet4p = dynamic_cast<StBET4pMaker*>(fourPMaker);
00087   if (bet4p) {
00088     jets.setDylanPoints( bet4p->nDylanPoints() );
00089     jets.setSumEmcE( bet4p->sumEmcEt() );
00090   }
00091 }
00092 
00093 void StjeDefaultJetTreeWriter::fillJetTree(int iAnalyzer, int iVertex)
00094 {
00095   if (iVertex) return;          // Only use first vertex
00096   StFourPMaker* fourPMaker = _analyzerCtlList[iAnalyzer]._fourPMaker;
00097   std::list<StProtoJet>* protoJetList = _analyzerCtlList[iAnalyzer]._protoJetList;
00098   fillJetTreeForOneJetFindingAlgorithm(*_analyzerCtlList[iAnalyzer]._jets, protoJetList, fourPMaker);
00099 }
00100 
00101 void StjeDefaultJetTreeWriter::fillJetTreeForOneJetFindingAlgorithm(StJets& jets, std::list<StProtoJet>* protoJetList, StFourPMaker* fourPMaker)
00102 {
00103   for(list<StProtoJet>::iterator it = protoJetList->begin(); it != protoJetList->end(); ++it) {
00104     fillJet(jets, *it);
00105   }
00106 }
00107 
00108 void StjeDefaultJetTreeWriter::fillJet(StJets &stjets, StProtoJet& pj)
00109 {
00110   StJet aJet(StJet(pj.e(), pj.px(), pj.py(), pj.pz(), 0, 0));
00111   stjets.addJet(aJet);
00112   TClonesArray* jets = stjets.jets();
00113   StJet* jet = (StJet*)jets->Last();
00114   jet->zVertex = _uDstMaker.muDst()->event()->primaryVertexPosition().z();
00115 
00116   const StProtoJet::FourVecList &particleList = pj.list();
00117   for(StProtoJet::FourVecList::const_iterator it2 = particleList.begin(); it2 != particleList.end(); ++it2)  {
00118     const StMuTrackFourVec *particle = dynamic_cast<const StMuTrackFourVec*>(*it2);
00119     if (!particle) {
00120       cout <<"StJets::addProtoJet(). ERROR:\tcast to StMuTrackFourVecFailed.  no action"<<endl;
00121       return;
00122     }
00123     int muTrackIndex = particle->getIndex();
00124     if (muTrackIndex <0) {
00125       cout <<"Error, muTrackIndex<0. abort()"<<endl;
00126       abort();
00127     }
00128 
00129     StDetectorId detectorId;
00130     int mDetId = particle->detectorId();
00131     if (mDetId==kTpcIdentifier)
00132       detectorId = kTpcId;
00133     else if (mDetId==kBarrelEmcTowerIdentifier)
00134       detectorId = kBarrelEmcTowerId;
00135     else if (mDetId==kEndcapEmcTowerIdentifier)
00136       detectorId = kEndcapEmcTowerId;
00137     else
00138       detectorId = kUnknownId;
00139 
00140     // Add tracks
00141     if (StMuTrackEmu* track = particle->track()) {
00142       TrackToJetIndex t2j( jets->GetLast(), muTrackIndex, detectorId, jet );
00143 
00144       t2j.SetPxPyPzE(particle->px(), particle->py(), particle->pz(), particle->e() );
00145       t2j.setTrackId( track->id() );
00146       t2j.setFlag( track->flag() );
00147       t2j.setCharge( track->charge() );
00148       t2j.setNhits( track->nHits() );
00149       t2j.setNhitsPoss( track->nHitsPoss() );
00150       t2j.setNhitsDedx( track->nHitsDedx() );
00151       t2j.setNhitsFit( track->nHitsFit() );
00152       t2j.setNsigmaPion( track->nSigmaPion() );
00153       t2j.setNsigmaElectron( track->nSigmaElectron() );
00154       t2j.setNsigmaKaon( track->nSigmaKaon() );
00155       t2j.setNsigmaProton( track->nSigmaProton() );
00156       t2j.setTdca ( track->Tdca() );
00157       t2j.setTdcaz ( track->dcaZ() );
00158       t2j.setTdcaxy ( track->dcaD() );
00159       t2j.setetaext ( track->etaext() );
00160       t2j.setphiext ( track->phiext() );
00161       t2j.setdEdx ( track->dEdx() );
00162 
00163       stjets.addTrackToIndex(t2j); // for backward compatibility
00164       jet->addTrack((TrackToJetIndex*)stjets.tracks()->Last());
00165     }
00166 
00167     // Add towers
00168     if (StMuTowerEmu* tower = particle->tower()) {
00169       TowerToJetIndex t2j(jets->GetLast());
00170 
00171       t2j.SetPxPyPzE(particle->px(), particle->py(), particle->pz(), particle->e());
00172       t2j.setTowerId(tower->id());
00173       t2j.setDetectorId(tower->detectorId());
00174       t2j.setAdc(tower->adc());
00175       t2j.setPedestal(tower->pedestal());
00176       t2j.setRms(tower->rms());
00177       t2j.setStatus(tower->status());
00178 
00179       stjets.addTowerToIndex(t2j); // for backward compatibility
00180       jet->addTower((TowerToJetIndex*)stjets.towers()->Last());
00181     }
00182 
00183     // Add Pythia particles
00184     if (StMcTrackEmu* mctrack = particle->mctrack()) {
00185       TrackToJetIndex t2j( jets->GetLast(), muTrackIndex, detectorId, jet );
00186 
00187       t2j.SetPtEtaPhiE(mctrack->pt(),mctrack->eta(),mctrack->phi(),mctrack->e());
00188       t2j.setId(mctrack->id());
00189       t2j.setPdg(mctrack->pdg());
00190       t2j.setStatus(mctrack->status());
00191 
00192       stjets.addTrackToIndex(t2j); // for backward compatibility
00193       jet->addTrack((TrackToJetIndex*)stjets.tracks()->Last());
00194     }
00195 
00196     if (mDetId==kTpcIdentifier) {
00197       jet->nTracks++;
00198       jet->tpcEtSum += particle->eT();
00199     }
00200     else if (mDetId==kBarrelEmcTowerIdentifier) {
00201       jet->nBtowers++;
00202       jet->btowEtSum += particle->eT();
00203     }
00204     else if (mDetId==kEndcapEmcTowerIdentifier) {
00205       jet->nEtowers++;
00206       jet->etowEtSum += particle->eT();
00207     }
00208   }
00209 }