StRoot: StSpinPool/StChargedPionAnalysisMaker/StChargedPionHelpers.cxx Source File

00001 // $Id: StChargedPionHelpers.cxx,v 1.7 2010/07/16 20:24:08 rfatemi Exp $
00002 
00003 #include "TLorentzVector.h"
00004 
00005 #include "StChargedPionHelpers.h"
00006 
00007 #include "StChargedPionBaseEv.h"
00008 #include "StChargedPionEvent.h"
00009 #include "StChargedPionJet.h"
00010 #include "StChargedPionMcEvent.h"
00011 #include "StChargedPionTrack.h"
00012 #include "StChargedPionVertex.h"
00013 
00014 #include "StSpinPool/StJetSkimEvent/StJetSkimEvent.h"
00015 #include "StSpinPool/StJets/StJets.h"
00016 #include "StSpinPool/StJets/StJet.h"
00017 #include "StMiniMcEvent/StMiniMcEvent.h"
00018 #include "StMuDSTMaker/COMMON/StMuDst.h"
00019 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00020 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
00021 #include "StSpinPool/StMCAsymMaker/StPythiaEvent.h"
00022 
00023 ClassImp(StChargedPionHelpers)
00024 
00025 StChargedPionHelpers::~StChargedPionHelpers() { }
00026 
00027 void StChargedPionHelpers::
00028 translateEvent(StJetSkimEvent *skimEvent, StChargedPionBaseEv *ev) {
00029     ev->setRunId(skimEvent->runId());
00030     ev->setEventId(skimEvent->eventId());
00031     ev->setBbcTimeBin(skimEvent->bbcTimeBin());
00032     char *baseName = strrchr(skimEvent->mudstFileName().GetString().Data(), '/');
00033     baseName = baseName + 1;
00034     ev->setMuDstName( baseName );
00035     
00036     StChargedPionEvent *data = dynamic_cast<StChargedPionEvent*>(ev);
00037     if(data) {
00038         data->setBx7(skimEvent->bx7());
00039 
00040         data->setSpinBit(skimEvent->spin4usingBx48());
00041         data->setPolValid(skimEvent->isValid());
00042         data->setPolLong(skimEvent->isPolLong());
00043         data->setPolTrans(skimEvent->isPolTrans());
00044         data->setBxingMasked(skimEvent->isMaskedUsingBx48());
00045         data->setBxingOffset(skimEvent->offsetBx48minusBX7());
00046 
00047         for(int i=0; i<skimEvent->trigHeaders()->GetEntriesFast(); i++) {
00048             StJetSkimTrigHeader *h = static_cast<StJetSkimTrigHeader*>(skimEvent->trigHeaders()->At(i));
00049             data->setPrescale(h->trigId, h->prescale);
00050         }
00051 
00052         for(int i=0; i<skimEvent->triggers()->GetEntriesFast(); i++) {
00053             StJetSkimTrig *t = static_cast<StJetSkimTrig*>(skimEvent->triggers()->At(i));
00054             if(t->didFire() > 0) data->addTrigger(t->trigId());
00055             if(t->shouldFire() > 0) data->addSimuTrigger(t->trigId());
00056 
00057             map<int,int> *m;
00058             m = &(t->towersAboveThreshold(0));
00059             for(map<int,int>::const_iterator it=m->begin(); it!=m->end(); it++) {
00060                 data->addHighTower(it->first, it->second);
00061             }
00062             delete m;
00063 
00064             m = &(t->triggerPatchesAboveThreshold(0));
00065             for(map<int,int>::const_iterator it=m->begin(); it!=m->end(); it++) {
00066                 data->addTriggerPatch(it->first, it->second);
00067             }
00068             delete m;
00069 
00070             m = &(t->jetPatchesAboveThreshold(0));
00071             for(map<int,int>::const_iterator it=m->begin(); it!=m->end(); it++) {
00072                 data->addJetPatch(it->first, it->second);
00073             }
00074             delete m;
00075         }        
00076     }
00077     
00078     if(skimEvent->eBbc() > 0 && skimEvent->wBbc() > 0) {
00079         if(ev->runId() > 6000000 && ev->runId() < 7000000) {
00080             ev->addSimuTrigger(96011);
00081         }
00082         else if(ev->runId() > 7000000 && ev->runId() < 8000000) {
00083             ev->addSimuTrigger(117001);
00084         }
00085     }
00086     
00087     StChargedPionVertex v;
00088     for(int i=0; i<skimEvent->vertices()->GetEntriesFast(); i++) {
00089         StJetSkimVert *skv = static_cast<StJetSkimVert*>(skimEvent->vertices()->At(i));
00090         v.SetX( skv->position()[0] );
00091         v.SetY( skv->position()[1] );
00092         v.SetZ( skv->position()[2] );
00093         v.setRanking( skv->ranking() );
00094         ev->addVertex(&v);
00095     }
00096 }
00097 
00098 void StChargedPionHelpers::
00099 translateJet(StJet* oldJet, vector<TLorentzVector*> particles, StChargedPionJet* jet) {
00100     *static_cast<TLorentzVector*>(jet) = *static_cast<TLorentzVector*>(oldJet);
00101     
00102     jet->setCharge(oldJet->charge);
00103     jet->setNTpcTracks(oldJet->nTracks);
00104     jet->setNBarrelTowers(oldJet->nBtowers);
00105     jet->setNEndcapTowers(oldJet->nEtowers);
00106     jet->setTpcEtSum(oldJet->tpcEtSum);
00107     jet->setBarrelEtSum(oldJet->btowEtSum);
00108     jet->setEndcapEtSum(oldJet->etowEtSum);
00109     jet->setVertexZ(oldJet->zVertex);
00110     
00111     for(unsigned i=0; i<oldJet->geomTriggers().size(); i++) {
00112         jet->addTrigger(oldJet->geomTriggers()[i]);
00113     }
00114     
00115     StChargedPionJetParticle particle;
00116     StThreeVectorF dca;
00117     TowerToJetIndex *myTower;
00118     TrackToJetIndex *myTrack;
00119     for( unsigned i=0; i<particles.size(); i++ )
00120     {
00121         particle.SetPt ( particles[i]->Pt()  );
00122         particle.SetEta( particles[i]->Eta() );
00123         particle.SetPhi( particles[i]->Phi() );
00124         particle.SetE( particles[i]->E() );
00125 
00126         if ( (myTower = dynamic_cast<TowerToJetIndex *>(particles[i])) )
00127         {
00128             particle.setIndex( myTower->towerId() );
00129             particle.setDetectorId( StDetectorId(myTower->detectorId()) );
00130             particle.setCharge( 0 );
00131         }
00132         else if ( (myTrack = dynamic_cast<TrackToJetIndex *>(particles[i])) )
00133         {
00134             particle.setIndex( myTrack->trackIndex() );
00135             particle.setDetectorId( myTrack->detectorId() );
00136             particle.setCharge( myTrack->charge());
00137             particle.setNHits( myTrack->nHits() );
00138             particle.setNHitsFit( myTrack->nHitsFit() );
00139             particle.setNHitsPoss( myTrack->nHitsPoss() );
00140             particle.setNHitsDEdx( myTrack->nHitsDedx() );
00141             particle.setNSigmaPion( myTrack->nSigmaPion() );
00142             
00143             dca.setY( float(myTrack->Tdcaxy()) );
00144             dca.setZ( float(myTrack->Tdcaz()) );
00145             particle.setGlobalDca( dca );
00146         }
00147         jet->addParticle(&particle);
00148     }
00149 }
00150 
00151 void StChargedPionHelpers::
00152 translateJets(StJets *jets, StChargedPionBaseEv *ev) {
00153     StChargedPionJet myjet;
00154     StChargedPionVertex *v = ev->vertex(0);
00155     
00156     TClonesArray *oldJets = jets->jets();
00157     for(int i=0; i<oldJets->GetEntriesFast(); i++) {
00158         StJet *oldJet = static_cast<StJet*>(oldJets->At(i));
00159         translateJet(oldJet, jets->particles(i), &myjet);
00160         if(v) myjet.setVertexZ(v->z());
00161         ev->addJet(&myjet);
00162         myjet.Clear();
00163     }
00164 }
00165 
00166 void StChargedPionHelpers::
00167 translateJets(StJets *jets, StChargedPionMcEvent *ev) {
00168     StChargedPionMcJet myjet;
00169     StChargedPionVertex *v = ev->vertex(0);
00170     
00171     TClonesArray *oldJets = jets->jets();
00172     for(int i=0; i<oldJets->GetEntriesFast(); i++) {
00173         StJet *oldJet = static_cast<StJet*>(oldJets->At(i));
00174         translateJet(oldJet, jets->particles(i), &myjet);
00175         if(v) myjet.setVertexZ(v->z());
00176         ev->mcJets().push_back(myjet);
00177         myjet.Clear();
00178     }
00179 }
00180 
00181 void StChargedPionHelpers::
00182 translatePythia(const StPythiaEvent *py, StChargedPionMcEvent *ev) {
00183     //remove const-ness until we have const versions of pythia->particle()
00184     StPythiaEvent *pythia = const_cast<StPythiaEvent*>(py);
00185     
00186     ev->parton1() = StChargedPionLorentzVector(pythia->particle(0)->Px(),
00187                                                pythia->particle(0)->Py(),
00188                                                pythia->particle(0)->Pz(),
00189                                                pythia->particle(0)->Energy());
00190     ev->parton2() = StChargedPionLorentzVector(pythia->particle(1)->Px(),
00191                                                pythia->particle(1)->Py(),
00192                                                pythia->particle(1)->Pz(),
00193                                                pythia->particle(1)->Energy());
00194     ev->parton3() = StChargedPionLorentzVector(pythia->particle(2)->Px(),
00195                                                pythia->particle(2)->Py(),
00196                                                pythia->particle(2)->Pz(),
00197                                                pythia->particle(2)->Energy());
00198     ev->parton4() = StChargedPionLorentzVector(pythia->particle(3)->Px(),
00199                                                pythia->particle(3)->Py(),
00200                                                pythia->particle(3)->Pz(),
00201                                                pythia->particle(3)->Energy());
00202     ev->setProcessId(pythia->processId());
00203     ev->setX1(pythia->x1());
00204 }
00205 
00206 void StChargedPionHelpers::
00207 translateMinimc(const StMiniMcEvent *mini, StChargedPionMcEvent *ev) {
00208     ev->mcVertex() = StChargedPion3Vector(mini->mcVertexX(),
00209                                           mini->mcVertexY(),
00210                                           mini->mcVertexZ());
00211     
00212     TClonesArray* mcTracks      = mini->tracks(MC);
00213     TClonesArray* matchedPairs  = mini->tracks(MATCHED);
00214     TClonesArray* mergedPairs   = mini->tracks(MERGED);
00215     TClonesArray* splitPairs    = mini->tracks(SPLIT);
00216     TClonesArray* contamPairs   = mini->tracks(CONTAM);
00217     TClonesArray* ghostPairs    = mini->tracks(GHOST);
00218 
00219     StTinyMcTrack *t;
00220     StMiniMcPair *p;
00221     for(int i=0; i<mcTracks->GetEntriesFast(); i++) {
00222         t = static_cast<StTinyMcTrack*>(mcTracks->At(i));
00223         if(t->ptMc() > 1.7) ev->mcTracks().push_back(*t);
00224     }
00225     
00226     //note: the approach I used here should end up duplicating the tiny mc tracks
00227     //in each of these clones arrays.  Alternative would be to find the pointer to
00228     //the original tiny mc track in ev->mcTracks() that has the same ID.  That's
00229     //likely a CPU-intensive process (but who cares?)
00230     
00231     //recoKey wasn't filled for the 2005 simulations, so I have to fall back to
00232     //pT matching here
00233     
00234     for(int i=0; i<matchedPairs->GetEntriesFast(); i++) {
00235         p = static_cast<StMiniMcPair*>(matchedPairs->At(i));
00236         if(p->ptPr() > 2.0) {
00237             ev->matchedPairs().push_back(*p);
00238         }
00239     }
00240     
00241     for(int i=0; i<mergedPairs->GetEntriesFast(); i++) {
00242         p = static_cast<StMiniMcPair*>(mergedPairs->At(i));
00243         if(p->ptPr() > 2.0) {
00244             ev->mergedPairs().push_back(*p);
00245         }
00246     }
00247     
00248     for(int i=0; i<splitPairs->GetEntriesFast(); i++) {
00249         p = static_cast<StMiniMcPair*>(splitPairs->At(i));
00250         if(p->ptPr() > 2.0) {
00251             ev->splitPairs().push_back(*p);
00252         }
00253     }
00254     
00255     for(int i=0; i<contamPairs->GetEntriesFast(); i++) {
00256         p = static_cast<StMiniMcPair*>(contamPairs->At(i));
00257         if(p->ptPr() > 2.0) {
00258             ev->contamPairs().push_back(*p);
00259         }
00260     }
00261     
00262     for(int i=0; i<ghostPairs->GetEntriesFast(); i++) {
00263         p = static_cast<StMiniMcPair*>(ghostPairs->At(i));
00264         if(p->ptPr() > 2.0) {
00265             ev->ghostPairs().push_back(*p);
00266         }
00267     }
00268 }
00269 
00270 void StChargedPionHelpers::
00271 translateMuDst(StChargedPionEvent *ev) {
00272     translateMuDst(static_cast<StChargedPionBaseEv*>(ev));
00273     
00274     ev->setBx7( StMuDst::event()->l0Trigger().bunchCrossingId7bit(ev->runId()) );
00275     ev->setRunInfo( StMuDst::event()->runInfo() );
00276     
00277     // see StDaqLib/TRG/trgStructures2005.h => L2RESULTS_OFFSET_DIJET==14
00278     if(ev->runId() > 7000000) 
00279         ev->setL2Result( StMuDst::event()->L2Result().GetArray() + 14 );
00280 }
00281 
00282 void StChargedPionHelpers::
00283 translateMuDst(StChargedPionMcEvent *ev) {
00284     translateMuDst(static_cast<StChargedPionBaseEv*>(ev));
00285 }
00286 
00287 void StChargedPionHelpers::
00288 translateMuDst(StChargedPionBaseEv *ev) {
00289     //basic event info
00290     ev->setRunId( StMuDst::event()->runNumber() );
00291     ev->setEventId( StMuDst::event()->eventNumber() );
00292     ev->setBbcTimeBin( StMuDst::event()->bbcTriggerDetector().onlineTimeDifference() );
00293     
00294     //vertices
00295     StChargedPionVertex vertex;
00296     unsigned int nVertices = StMuDst::numberOfPrimaryVertices();
00297     for(unsigned int i=0; i<nVertices; ++i){
00298         assert(StMuDst::primaryVertex(i));
00299         StMuPrimaryVertex* muVert = StMuDst::primaryVertex(i);
00300         vertex.SetX( muVert->position().x() );
00301         vertex.SetY( muVert->position().y() );
00302         vertex.SetZ( muVert->position().z() );
00303         vertex.setRanking( muVert->ranking() );
00304         ev->addVertex(&vertex);
00305     }
00306     
00307     //now for the tracks
00308     for(unsigned int vertex_index=0; vertex_index<nVertices; vertex_index++){
00309         StMuDst::setVertexIndex(vertex_index);
00310         TObjArray* primaryTracks = StMuDst::primaryTracks();
00311         const StMuTrack* track;
00312         const StMuTrack* global;
00313         int nentries = StMuDst::numberOfPrimaryTracks();
00314         assert(nentries==primaryTracks->GetEntries());
00315         for(int i=0; i<nentries; i++){
00316             track = StMuDst::primaryTracks(i);
00317             global = track->globalTrack();
00318             
00319             //sanity checks
00320             switch(track->bad()){
00321                 case(0):
00322                     break;
00323                 case(10):
00324                     // mBadTracks->Fill(kFlagged);
00325                 case(20):
00326                     // mBadTracks->Fill(kBadHelix);
00327                 case(30):
00328                     // mBadTracks->Fill(kBadOuterHelix);
00329                     continue;
00330             }
00331             
00332             if(!global){
00333                 // mBadTracks->Fill(kMissingGlobal);
00334                 continue;
00335             }
00336             
00337             //cuts
00338             if(track->pt() < 2.)                continue;
00339             if(TMath::Abs(track->eta()) > 2.1)  continue;
00340             
00341             StChargedPionTrack cpTrack;
00342             translateTrack(track, &cpTrack);
00343             cpTrack.setB(StMuDst::event()->eventSummary().magneticField()*kilogauss);
00344             cpTrack.setVertex(StMuDst::primaryVertex()->position());
00345             
00346             ev->addTrack(&cpTrack);
00347         }
00348     }
00349 }
00350 
00351 void StChargedPionHelpers::
00352 translateTrack(const StMuTrack *mu, StChargedPionTrack *cp) {
00353     cp->setId(mu->id());
00354     cp->setFlag(mu->flag());
00355     
00356     cp->setVertexIndex(mu->vertexIndex());
00357     
00358     cp->setNHits(mu->nHits());
00359     cp->setNHitsPoss(mu->nHitsPoss());
00360     cp->setNHitsDedx(mu->nHitsDedx());
00361     cp->setNHitsFit(mu->nHitsFit());
00362     
00363     cp->setNSigmaElectron(mu->nSigmaElectron());
00364     cp->setNSigmaPion(mu->nSigmaPion());
00365     cp->setNSigmaKaon(mu->nSigmaKaon());
00366     cp->setNSigmaProton(mu->nSigmaProton());
00367     cp->setDedx(mu->dEdx());
00368     
00369     cp->setChi2(mu->chi2());
00370     cp->setChi2prob(mu->chi2prob());
00371     
00372     cp->setPtEtaPhi(mu->pt(),mu->eta(),mu->phi());
00373     
00374     cp->setCharge(mu->charge());      
00375     
00376     cp->setGlobalLastPoint(mu->globalTrack()->lastPoint());
00377     cp->setGlobalHelix(mu->globalTrack()->helix());
00378 }
00379 
00380 /*****************************************************************************
00381  * $Log: StChargedPionHelpers.cxx,v $
00382  * Revision 1.7  2010/07/16 20:24:08  rfatemi
00383  * Changes in code to preserve backward compatibility with changes in jet code.  Affects info about particles in jets.
00384  *
00385  * Revision 1.5  2009/12/02 21:17:31  fine
00386  * Fix StMuTrack interface
00387  *
00388  * Revision 1.4  2009/09/25 14:26:37  fine
00389  * fix compilation error on SL5/64-bit machine
00390  *
00391  * Revision 1.3  2009/04/02 18:25:42  kocolosk
00392  * fixed paths to jet codes
00393  *
00394  * Revision 1.2  2008/12/29 15:58:28  kocolosk
00395  * removed commented code and added Id and Log as needed
00396  *
00397  * Revision 1.1  2008/07/17 17:06:30  kocolosk
00398  * big-bang integration StChargedPionMcEvent framework
00399  *
00400  *****************************************************************************/