StRoot: StSpinPool/StChargedPionAnalysisMaker/StChargedPionMaker.cxx Source File

00001 // $Id: StChargedPionMaker.cxx,v 1.22 2010/04/25 16:07:57 pibero Exp $
00002 
00003 #include "StChargedPionMaker.h"
00004 
00005 //ROOT headers
00006 #include "TFile.h"
00007 #include "TTree.h"
00008 #include "TClonesArray.h"
00009 #include "TChain.h"
00010 
00011 //StarClassLibrary
00012 #include "SystemOfUnits.h"
00013 
00014 //StMuDstMaker
00015 #include "StMuDSTMaker/COMMON/StMuDst.h"
00016 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00017 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00018 #include "StMuDSTMaker/COMMON/StMuDst.h"
00019 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00020 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
00021 
00022 //logger
00023 #include "StMessMgr.h"
00024 
00025 //my headers
00026 #include "StChargedPionEvent.h"
00027 #include "StChargedPionVertex.h"
00028 #include "StChargedPionTrack.h"
00029 #include "StChargedPionJet.h"
00030 #include "StChargedPionJetParticle.h"
00031 #include "StChargedPionMcEvent.h"
00032 #include "StChargedPionHelpers.h"
00033 #include "StChargedPionTypes.h"
00034 
00035 //StDetectorDbMaker
00036 #include "StDetectorDbMaker/StDetectorDbTriggerID.h"
00037 
00038 //StJetMaker
00039 #include "StSpinPool/StJetSkimEvent/StJetSkimEvent.h"
00040 #include "StSpinPool/StJets/StJets.h"
00041 #include "StSpinPool/StJets/StJet.h"
00042 #include "StJetMaker/StJetMaker.h"
00043 
00044 //StSpinDbMaker
00045 #include "StSpinPool/StSpinDbMaker/StSpinDbMaker.h"
00046 
00047 //StEmcTriggerMaker
00048 #include "StEmcTriggerMaker/StEmcTriggerMaker.h"
00049 
00050 //StTriggerUtilities
00051 #include "StTriggerUtilities/StTriggerSimuMaker.h"
00052 #include "StTriggerUtilities/StTriggerSimuResult.h"
00053 
00054 //StMiniMcEvent
00055 #include "StMiniMcEvent/StMiniMcEvent.h"
00056 #include "StMiniMcEvent/StTinyMcTrack.h"
00057 #include "StMiniMcEvent/StMiniMcPair.h"
00058 
00059 #include "tables/St_g2t_event_Table.h"
00060 #include "tables/St_particle_Table.h"
00061 #include "tables/St_g2t_pythia_Table.h"
00062 
00063 //StMcEvent
00064 #include "StMcEvent/StMcEvent.hh"
00065 
00066 ClassImp(StChargedPionMaker)
00067 
00068 StChargedPionMaker::StChargedPionMaker(const char *name, const char *outputfile) : 
00069     StMaker(name) 
00070 {
00071     LOG_INFO << "calling constructor" << endm;
00072     
00073     mFile = new TFile(outputfile,"RECREATE");
00074 
00075     mBadTracks = new TH1D("badTracks","tracks failing quality cuts",4,0.5,4.5);
00076 
00077     time_t rawtime;
00078     time(&rawtime);
00079     char title[200];
00080     sprintf(title, "created %s", ctime(&rawtime));
00081     mTree = new TTree("tree",title);
00082     
00083     bool isMonteCarlo = GetMakerInheritsFrom("StEmcSimulatorMaker") ? true:false;
00084     LOG_INFO << "Are we running on Monte Carlo? " << isMonteCarlo << endm;
00085     
00086     if(isMonteCarlo) {
00087         mEvent = new StChargedPionMcEvent();
00088         mTree->Branch("event", "StChargedPionMcEvent", &mEvent);
00089     }
00090     else {
00091         mEvent = new StChargedPionEvent();
00092         mTree->Branch("event", "StChargedPionEvent", &mEvent);
00093     }
00094     
00095     Long64_t autosave = 1000000000; //1GB
00096     mTree->SetAutoSave(autosave);
00097     mTree->SetMaxTreeSize(autosave);
00098     
00099     mMuDstMk      = NULL;
00100     mSpDbMk       = NULL;
00101     mEmcTrgMk     = NULL;
00102     mJetMk        = NULL;
00103     mTrgSimuMk   = NULL;
00104     
00105     mJetFile = NULL;
00106     mJetTree = NULL;
00107     mJets = new StJets();
00108     mPyJets = new StJets();
00109     
00110     mMiniMcFile = NULL;
00111     mMiniMcTree = NULL;
00112     mMiniMcEvent = new StMiniMcEvent();
00113     
00114     mTriggers.clear();
00115     
00116     LOG_INFO << "finished constructor" << endm; 
00117 }
00118 
00119 StChargedPionMaker::~StChargedPionMaker() {
00120     LOG_DEBUG << "calling destructor" << endm;
00121     
00122     if(mEvent) delete mEvent;
00123     
00124     if(mJetFile) mJetFile->Close();
00125     delete mJets;
00126     delete mPyJets;
00127     
00128     if(mMiniMcFile) mMiniMcFile->Close();
00129     delete mMiniMcEvent;
00130     
00131     LOG_DEBUG << "finished destructor" << endm;
00132 }
00133 
00134 void StChargedPionMaker::Clear(const char*) {
00135     StChargedPionEvent   *data = dynamic_cast<StChargedPionEvent*>(mEvent);
00136     StChargedPionMcEvent *simu = dynamic_cast<StChargedPionMcEvent*>(mEvent);
00137     if(data) data->Clear();
00138     if(simu) simu->Clear();
00139     StMaker::Clear();
00140 }
00141 
00142 Int_t StChargedPionMaker::Init() {
00143     mMuDstMk    = dynamic_cast<StMuDstMaker*>(
00144                         GetMakerInheritsFrom("StMuDstMaker"));
00145     mSpDbMk     = dynamic_cast<StSpinDbMaker*>(
00146                         GetMakerInheritsFrom("StSpinDbMaker"));
00147     mEmcTrgMk   = dynamic_cast<StEmcTriggerMaker*>(
00148                         GetMakerInheritsFrom("StEmcTriggerMaker")); 
00149     mJetMk      = dynamic_cast<StJetMaker*>(
00150                         GetMakerInheritsFrom("StJetMaker"));
00151     mTrgSimuMk  = dynamic_cast<StTriggerSimuMaker*>(
00152                         GetMakerInheritsFrom("StTriggerSimuMaker"));
00153     
00154     // do the setup for mPyJets here, since we only run on one file at a time
00155     if( dynamic_cast<StChargedPionMcEvent*>(mEvent) ) {
00156         if(mJetMk) {
00157             mPyJets = mJetMk->getStJets("PythiaConeJets");
00158             LOG_INFO << "loaded mPyJets at " << mPyJets << endm;
00159         }
00160     }
00161     
00162     this->Clear();
00163     
00164     LOG_INFO << "init OK" << endm;
00165     return StMaker::Init();
00166 }
00167 
00168 Int_t StChargedPionMaker::InitRun(int runnumber) {
00169     if(mJetMk) {
00170         LOG_INFO << "found StJetMaker in the chain" << endm;
00171         if(StJets* stjets = mJetMk->getStJets("ConeJets")) {
00172             mJets = stjets;
00173             LOG_INFO << "found Jets in Run 5 branch " << mJets << endm;
00174         }
00175         else if(StJets* stjets = mJetMk->getStJets("ConeJets12")) {
00176             mJets = stjets;
00177             LOG_INFO << "found Jets in Run 6 branch " << mJets << endm;
00178         }
00179         else {
00180             mJets = mJetMk->getStJets("ConeJets12_0.7");
00181             LOG_INFO << "found ConeJets12_0.7 branch " << mJets << endm;
00182         }
00183     }
00184     else {
00185         LOG_INFO << "trying to get the jets off disk" << endm;
00186         std::ostringstream os;
00187         if(runnumber < 7000000) {
00188             os << "/star/institutions/mit/common/run5/jets/jets_" << runnumber << ".tree.root";
00189         }
00190         else {
00191             // best confirm this one with Murad first
00192             os << "/star/institutions/mit/common/run6/jets/jets_" << runnumber << ".tree.root";
00193         }
00194         
00195         if(mJetFile) mJetFile->Close();
00196         mJetTree = NULL;
00197         mJetFile = TFile::Open(os.str().c_str());
00198         if(mJetFile) mJetTree = (TTree*) mJetFile->Get("jet");
00199         if(mJetTree) {
00200             if(runnumber < 7000000) {
00201                 mJetTree->SetBranchAddress("ConeJets", &mJets);            
00202             }
00203             else {
00204                 mJetTree->SetBranchAddress("ConeJets12", &mJets);            
00205             }
00206             mJetTree->BuildIndex("mRunId","mEventId");
00207         }
00208     }
00209     
00210     return StMaker::InitRun(runnumber);
00211 }
00212 
00213 Int_t StChargedPionMaker::Make()
00214 {
00215     StChargedPionEvent   *data = dynamic_cast<StChargedPionEvent*>(mEvent);
00216     StChargedPionMcEvent *simu = dynamic_cast<StChargedPionMcEvent*>(mEvent);
00217     
00218     //have we changed files?
00219     TString inputFile(mMuDstMk->chain()->GetFile()->GetName());
00220     if(mCurrentFile != inputFile){
00221         mCurrentFile = inputFile;
00222         char *baseName = strrchr(mCurrentFile.Data(), '/');
00223         mEvent->setMuDstName( baseName );
00224         
00225         if(simu) {
00226             TString minimcFile = inputFile.ReplaceAll("MuDst","minimc");
00227             mMiniMcFile = TFile::Open(minimcFile);
00228             if(mMiniMcFile) {
00229                 LOG_INFO << "opened minimc file at " << minimcFile << endm;
00230                 mMiniMcTree = dynamic_cast<TTree*> (mMiniMcFile->Get("StMiniMcTree"));
00231                 mMiniMcTree->BuildIndex("mEventId");
00232                 mMiniMcTree->SetBranchAddress("StMiniMcEvent", &mMiniMcEvent);
00233             }
00234             else {
00235                 LOG_WARN << "problem opening minimc at " << minimcFile << endm;
00236             }
00237         }
00238     }
00239     
00240     if(data) {
00241         StChargedPionHelpers::translateMuDst(data);
00242         
00243         // triggers, prescales
00244         mTriggers.clear();
00245         map<int,float> m = StDetectorDbTriggerID::instance()->getTotalPrescales();
00246         for (map<int,float>::iterator it=m.begin(); it!=m.end(); ++it) {
00247             int trigId = it->first;
00248             data->setPrescale(trigId, it->second);
00249                         
00250             if( StMuDst::event()->triggerIdCollection().nominal().isTrigger(trigId) ) {
00251                 data->addTrigger(trigId);
00252             }
00253             
00254             mTriggers.push_back(trigId);
00255         }
00256         
00257         makeTriggerSimu(data);
00258         
00259         //spin DB
00260         int bx48 =  StMuDst::event()->l0Trigger().bunchCrossingId();
00261         data->setSpinBit( mSpDbMk->spin4usingBX48(bx48) );
00262         data->setPolValid( mSpDbMk->isValid() );
00263         data->setPolLong( mSpDbMk->isPolDirLong() );
00264         data->setPolTrans( mSpDbMk->isPolDirTrans() );
00265         data->setBxingMasked( mSpDbMk->isMaskedUsingBX48(bx48) );
00266         data->setBxingOffset( mSpDbMk->offsetBX48minusBX7(bx48, data->bx7()) );
00267     }
00268     
00269     if(simu) {
00270         StChargedPionHelpers::translateMuDst(simu);
00271         
00272         makeTriggerSimu(simu);
00273         
00274         StMcEvent *mcEvent = static_cast<StMcEvent*>( GetDataSet("StMcEvent") );
00275         simu->setProcessId( mcEvent->subProcessId() );
00276         
00277         TDataSet *Event = GetDataSet("geant");
00278         TDataSetIter geantDstI(Event);
00279 
00280         St_g2t_pythia* PyPtr = (St_g2t_pythia *) geantDstI("g2t_pythia");
00281         g2t_pythia_st *g2t_pythia = PyPtr->GetTable();
00282         simu->setHardP( g2t_pythia->hard_p );
00283         simu->setX1( g2t_pythia->bjor_1 );
00284         
00285         St_particle *particleTabPtr = (St_particle *) geantDstI("particle");
00286         particle_st* particleTable = particleTabPtr->GetTable();
00287         
00288         simu->isr1().SetXYZT( 
00289             particleTable[2].phep[0],
00290             particleTable[2].phep[1],
00291             particleTable[2].phep[2],
00292             particleTable[2].phep[3]
00293             );
00294         
00295         simu->isr2().SetXYZT(
00296             particleTable[3].phep[0],
00297             particleTable[3].phep[1],
00298             particleTable[3].phep[2],
00299             particleTable[3].phep[3]
00300             );
00301         
00302         simu->parton1().SetXYZT(
00303             particleTable[4].phep[0],
00304             particleTable[4].phep[1],
00305             particleTable[4].phep[2],
00306             particleTable[4].phep[3]
00307             );
00308         
00309         simu->parton2().SetXYZT(
00310             particleTable[5].phep[0],
00311             particleTable[5].phep[1],
00312             particleTable[5].phep[2],
00313             particleTable[5].phep[3]
00314             );
00315 
00316         simu->parton3().SetXYZT(
00317             particleTable[6].phep[0],
00318             particleTable[6].phep[1],
00319             particleTable[6].phep[2],
00320             particleTable[6].phep[3]
00321             );
00322 
00323         simu->parton4().SetXYZT(
00324             particleTable[7].phep[0],
00325             particleTable[7].phep[1],
00326             particleTable[7].phep[2],
00327             particleTable[7].phep[3]
00328             );
00329 
00330         simu->setFlavor(1, particleTable[4].idhep);
00331         simu->setFlavor(2, particleTable[5].idhep);
00332         simu->setFlavor(3, particleTable[6].idhep);
00333         simu->setFlavor(4, particleTable[7].idhep);
00334         
00335         // save all final-state particles with E > 2.0 GeV
00336         for(int i = 8; i < particleTabPtr->GetNRows(); ++i)
00337         {
00338             if(particleTable[i].isthep == 1 && particleTable[i].phep[3] > 2.0)
00339             {
00340                 StChargedPionPythiaRow r;
00341                 r.id = particleTable[i].idhep;
00342                 r.vec = StChargedPionLorentzVector(particleTable[i].phep[0],
00343                                                    particleTable[i].phep[1],
00344                                                    particleTable[i].phep[2],
00345                                                    particleTable[i].phep[3]);
00346                 simu->pythiaRecord().push_back(r);
00347             }
00348         }
00349 
00350         int bytesRead = mMiniMcTree->GetEntryWithIndex(simu->eventId());
00351         if(bytesRead > 0) {
00352             StChargedPionHelpers::translateMinimc(mMiniMcEvent, simu);
00353         }
00354 
00355         // still need to set the the pythia jets
00356         if(mJetMk) {
00357             StChargedPionHelpers::translateJets(mPyJets, simu);
00358         }
00359         else if(mJetTree) {
00360             int ok = mJetTree->GetEntryWithIndex(simu->runId(), simu->eventId());
00361             if(ok > 0) StChargedPionHelpers::translateJets(mJets, simu);
00362         }
00363     }
00364 
00365     //and the jets
00366     if(mJetMk) {
00367         StChargedPionHelpers::translateJets(mJets, mEvent);
00368     }
00369     else if(mJetTree) {
00370         int ok = mJetTree->GetEntryWithIndex(mEvent->runId(), mEvent->eventId());
00371         if(ok > 0) StChargedPionHelpers::translateJets(mJets, mEvent);
00372     }
00373     
00374     mTree->Fill();
00375     
00376     return StMaker::Make();
00377 }
00378 
00379 Int_t StChargedPionMaker::Finish() {
00380     mFile->cd();
00381     mTree->Write();
00382     mBadTracks->Write();
00383     mFile->Close();
00384     LOG_INFO << "finished OK"<<endm;
00385     return StMaker::Finish();
00386 }
00387 
00388 void StChargedPionMaker::makeTriggerSimu(StChargedPionBaseEv *ev) {
00389     if(mTrgSimuMk) {
00390         for (unsigned i=0; i<mTriggers.size(); ++i) {
00391             int trigId = mTriggers[i];
00392             
00393             if( mTrgSimuMk->isTrigger(trigId) ) {
00394                 ev->addSimuTrigger(trigId);
00395             }
00396             
00397             const StTriggerSimuResult result = mTrgSimuMk->detailedResult(trigId);
00398             for(unsigned i=0; i<result.highTowerIds().size(); i++) {
00399                 int tid = result.highTowerIds().at(i);
00400                 ev->addHighTower(tid, result.highTowerAdc(tid));
00401             }
00402             for(unsigned i=0; i<result.triggerPatchIds().size(); i++) {
00403                 int tid = result.triggerPatchIds().at(i);
00404                 ev->addTriggerPatch(tid, result.triggerPatchAdc(tid));
00405             }
00406             for(unsigned i=0; i<result.jetPatchIds().size(); i++) {
00407                 int tid = result.jetPatchIds().at(i);
00408                 ev->addJetPatch(tid, result.jetPatchAdc(tid));
00409             }
00410             ev->setL2Result(result.l2Result(kJet), true);
00411         }
00412     }
00413     else if(mEmcTrgMk) {
00414         //minbias simu trigger
00415         int Npmt=StMuDst::event()->bbcTriggerDetector().numberOfPMTs();
00416         bool eastBBC(false), westBBC(false);
00417         for (int pmt=0;pmt<Npmt;pmt++){
00418             if(StMuDst::event()->bbcTriggerDetector().adc(pmt) > 5) {
00419                 if(pmt<16) eastBBC = true;
00420                 if(pmt>23 && pmt<40) westBBC = true;
00421             }
00422         }
00423         if(eastBBC && westBBC) {
00424             ev->addSimuTrigger(96011);  
00425             ev->addSimuTrigger(117011);
00426         }
00427         
00428         for (unsigned i=0; i<mTriggers.size(); ++i) {
00429             int trigId = mTriggers[i];
00430             
00431             if ( mEmcTrgMk->isTrigger(trigId) ) {
00432                 ev->addSimuTrigger(trigId);
00433             }
00434             
00435             map<int,int> m( mEmcTrgMk->barrelTowersAboveThreshold(trigId) );
00436             for(map<int,int>::const_iterator iter=m.begin(); iter!=m.end(); iter++) {
00437                 ev->addHighTower(iter->first, iter->second);
00438             }
00439             
00440             m = mEmcTrgMk->barrelTriggerPatchesAboveThreshold(trigId);
00441             for(map<int,int>::const_iterator iter=m.begin(); iter!=m.end(); iter++) {
00442                 ev->addTriggerPatch(iter->first, iter->second);
00443             }
00444             
00445             m = mEmcTrgMk->barrelJetPatchesAboveThreshold(trigId);
00446             for(map<int,int>::const_iterator iter=m.begin(); iter!=m.end(); iter++) {
00447                 ev->addJetPatch(iter->first, iter->second);
00448             }
00449         }
00450     }
00451 }
00452 
00453 /*****************************************************************************
00454  * $Log: StChargedPionMaker.cxx,v $
00455  * Revision 1.22  2010/04/25 16:07:57  pibero
00456  * Modified StChargedPionMaker.cxx to use the safer StJetMaker::getStJets()
00457  * instead of the outdated and no longer available StJetMaker::getJets().
00458  *
00459  * Revision 1.21  2010/04/25 15:18:49  pibero
00460  * Temporary fix to keep AutoBuild happy. Permanent solution to be posted soon.
00461  *
00462  * Revision 1.20  2009/06/22 14:44:50  kocolosk
00463  * support for new ConeJets12_0.7 branch name
00464  *
00465  * Revision 1.19  2009/04/02 18:25:42  kocolosk
00466  * fixed paths to jet codes
00467  *
00468  * Revision 1.18  2008/12/29 15:58:30  kocolosk
00469  * removed commented code and added Id and Log as needed
00470  *
00471  * Revision 1.17  2008/08/25 20:55:35  kocolosk
00472  * get correct ConeJets/ConeJets12 branch w/o using runnumber
00473  *
00474  * Revision 1.16  2008/07/17 17:06:30  kocolosk
00475  * big-bang integration StChargedPionMcEvent framework
00476  *
00477  *****************************************************************************/