StMcEvent.cc

#include "StThreeVectorF.hh"

#include "StMcEvent.hh"
#include "StMcTpcHitCollection.hh"
#include "StMcFtpcHitCollection.hh"
#include "StMcRichHitCollection.hh"
#include "StMcCtbHitCollection.hh"
#include "StMcSvtHitCollection.hh"
#include "StMcSsdHitCollection.hh"
#include "StMcEmcHitCollection.hh"
#include "StMcTofHitCollection.hh"
#include "StMcBTofHitCollection.hh"
#include "StMcMtdHitCollection.hh"
#include "StMcPxlHitCollection.hh"
#include "StMcIstHitCollection.hh"
#include "StMcFgtHitCollection.hh"
#include "StMcEtrHitCollection.hh"
#include "StMcContainers.hh" 
#include "StMcVertex.hh"
#include "StMcTrack.hh"
#include "StMcTpcHit.hh"
#include "StMcFtpcHit.hh"
#include "StMcRichHit.hh"
#include "StMcSvtHit.hh"
#include "StMcSsdHit.hh"
#include "StMcCalorimeterHit.hh"
#include "StMcTofHit.hh"
#include "StMcBTofHit.hh"
#include "StMcMtdHit.hh"
#include "StMcPxlHit.hh"
#include "StMcIstHit.hh"
#include "StMcFgtHit.hh"
#include "StMcEtrHit.hh"
#include "tables/St_g2t_event_Table.h"
#include "TDataSetIter.h"


TString StMcEvent::mCvsTag = "$Id: StMcEvent.cc,v 2.40 2016/05/17 19:41:50 perev Exp $";
static const char rcsid[] = "$Id: StMcEvent.cc,v 2.40 2016/05/17 19:41:50 perev Exp $";
ClassImp(StMcEvent);
//______________________________________________________________________________
void StMcEvent::initToZero()
{
    memset(mBeg,0,mEnd-mBeg+1);    
}
//______________________________________________________________________________
void StMcEvent::makeColls()
{
    // Create the collections
    
    mTpcHits  = new StMcTpcHitCollection();
    mSvtHits  = new StMcSvtHitCollection();
    mSsdHits  = new StMcSsdHitCollection();
    mFtpcHits = new StMcFtpcHitCollection();
    mRichHits = new StMcRichHitCollection();
    mCtbHits  = new StMcCtbHitCollection();
    mTofHits  = new StMcTofHitCollection();
    mBTofHits = new StMcBTofHitCollection();
    mMtdHits  = new StMcMtdHitCollection();
    mPxlHits  = new StMcPxlHitCollection();
    mIstHits  = new StMcIstHitCollection();
    mFgtHits  = new StMcFgtHitCollection();
    mEtrHits  = new StMcEtrHitCollection();
}

//______________________________________________________________________________
StMcEvent::StMcEvent() :TDataSet("StMcEvent")
{
    initToZero();
}

//______________________________________________________________________________
StMcEvent::StMcEvent(g2t_event_st* evTable)
    :TDataSet("StMcEvent")
{
     initToZero();
     mEventGeneratorEventLabel = evTable->eg_label;
     mEventNumber = evTable->n_event;
     mRunNumber = evTable->n_run;
     mType  = evTable->event_type;
     mZWest = evTable->n_part_prot_west;
     mNWest = evTable->n_part_neut_west;
     mZEast = evTable->n_part_prot_east;
     mNEast = evTable->n_part_neut_east;
     mEvGenFSTracks = evTable->n_track_eg_fs;
     mPrimaryTracks = evTable->n_track_prim;
     mSubProcessId = evTable->subprocess_id;
     mImpactParameter = evTable->b_impact;
     mPhiReactionPlane = evTable->phi_impact;
     mTriggerTimeOffset = evTable->time_offset;
     mNBinary = evTable->n_binary;
     mNWoundedEast = evTable->n_wounded_east;
     mNWoundedWest = evTable->n_wounded_west;
     mNJets = evTable->njets;
     makeColls();
}


//______________________________________________________________________________
StMcEvent::StMcEvent(const StMcEvent&) { /* noop */} // private

//______________________________________________________________________________
const StMcEvent& StMcEvent::operator=(const StMcEvent&) { return *this;} // private

//______________________________________________________________________________
StMcEvent::~StMcEvent()
{
  for (auto jk=mBegColl; jk<mEndColl;jk++) { 
    delete *jk; *jk=0;
  }
  for(auto jk=0;jk<(int)mTracks.size()  ;jk++) {delete mTracks[jk]  ;}    
  mTracks.clear();
  
  for(int jk=0;jk<(int)mVertices.size();jk++) {delete mVertices[jk];}
  mVertices.clear();
}

//______________________________________________________________________________
int StMcEvent::operator==(const StMcEvent& e) const
{
  return (e.mEventNumber == mEventNumber &&
          e.mRunNumber   == mRunNumber &&
          e.mType        == mType
          );  
}

//______________________________________________________________________________
int StMcEvent::operator!=(const StMcEvent& e) const
{
    return !(e == *this);   // invoke operator==()
}

//______________________________________________________________________________
ostream&  operator<<(ostream& os, const StMcEvent& e)
{
    os << "Label: " << e.eventGeneratorEventLabel()
       << "\tRun: " << e.runNumber()
       << "\tId: " << e.eventNumber() 
       << "\tType: " << e.type() << endl;
    os << "Participant Protons  East: " << e.zEast() << "\tWest: " << e.zWest() << endl;
    os << "Participant Neutrons East: " << e.nEast() << "\tWest: " << e.nWest() << endl;
    os << "# Ev. Gen. Fin. St. Track: " << e.eventGeneratorFinalStateTracks() 
       << "\tNumber of Primary Tracks : " << e.numberOfPrimaryTracks() << endl;
    os << "Subprocess Id    : " << e.subProcessId() << endl;
    os << "Impact Parameter : " << e.impactParameter() << endl;
    os << "Phi Reaction Pl. : " << e.phiReactionPlane() << endl;
    os << "Trig. Time Offset: " << e.triggerTimeOffset() << endl;
    os << "N Binary Coll.   : " << e.nBinary() << endl;
    os << "N Wounded East   : " << e.nWoundedEast() << endl;
    os << "N Wounded West   : " << e.nWoundedWest() << endl;
    os << "N Jets           : " << e.nJets() << endl;
    return os;
}


//______________________________________________________________________________
void StMcEvent::setEventGeneratorEventLabel(unsigned long val) { mEventGeneratorEventLabel = val; }

//______________________________________________________________________________
void StMcEvent::setEventNumber(unsigned long  val) { mEventNumber = val;  }

//______________________________________________________________________________
void StMcEvent::setRunNumber(unsigned long val) { mRunNumber = val; }                

//______________________________________________________________________________
void StMcEvent::setType(unsigned long val) { mType = val; }              

//______________________________________________________________________________
void StMcEvent::setZWest(unsigned long val) { mZWest = val; }              
//______________________________________________________________________________

void StMcEvent::setNWest(unsigned long val) { mNWest = val; }     

//______________________________________________________________________________
void StMcEvent::setZEast(unsigned long val) { mZEast = val; }              

//______________________________________________________________________________
void StMcEvent::setNEast(unsigned long val) { mNEast = val; }     

//______________________________________________________________________________
void StMcEvent::setEventGeneratorFinalStateTracks(unsigned long val){ mEvGenFSTracks = val; } 

//______________________________________________________________________________
void StMcEvent::setNumberOfPrimaryTracks(unsigned long val){ mPrimaryTracks = val; } 

//______________________________________________________________________________
void StMcEvent::setImpactParameter(float val) { mImpactParameter = val; }               

//______________________________________________________________________________
void StMcEvent::setPhiReactionPlane(float val) { mPhiReactionPlane = val; } 

//______________________________________________________________________________
void StMcEvent::setTriggerTimeOffset(float val) { mTriggerTimeOffset = val; }

//______________________________________________________________________________
void StMcEvent::setNBinary(unsigned long val) { mNBinary = val; }

//______________________________________________________________________________
void StMcEvent::setNWoundedEast(unsigned long val) { mNWoundedEast = val; }

//______________________________________________________________________________
void StMcEvent::setNWoundedWest(unsigned long val) { mNWoundedWest = val; }

//______________________________________________________________________________
void StMcEvent::setNJets(unsigned long val) { mNJets = val; }

//______________________________________________________________________________
void StMcEvent::setPrimaryVertex(StMcVertex* val) {  mPrimaryVertex = val; }

//______________________________________________________________________________
void StMcEvent::setTpcHitCollection(StMcTpcHitCollection* val)
{   
    if (mTpcHits!= val) delete mTpcHits;
    mTpcHits = val;
}               

//______________________________________________________________________________
void StMcEvent::setSvtHitCollection(StMcSvtHitCollection* val)
{
    if (mSvtHits!= val) delete mSvtHits;
    mSvtHits = val;
}               

//______________________________________________________________________________
void StMcEvent::setSsdHitCollection(StMcSsdHitCollection* val)
{
    if (mSsdHits!= val) delete mSsdHits;
    mSsdHits = val;
}               

//______________________________________________________________________________
void StMcEvent::setFtpcHitCollection(StMcFtpcHitCollection* val)
{
    if (mFtpcHits!= val) delete mFtpcHits;
    mFtpcHits = val;
}              

//______________________________________________________________________________
void StMcEvent::setRichHitCollection(StMcRichHitCollection* val)
{
    if (mRichHits!= val) delete mRichHits;
    mRichHits = val;
}              
//______________________________________________________________________________
void StMcEvent::setTofHitCollection(StMcTofHitCollection* val)
{
    if (mTofHits!= val) delete mTofHits;
    mTofHits = val;
}
//______________________________________________________________________________
void StMcEvent::setBTofHitCollection(StMcBTofHitCollection* val)
{
    if (mBTofHits!= val) delete mBTofHits;
    mBTofHits = val;
}
//______________________________________________________________________________
void StMcEvent::setMtdHitCollection(StMcMtdHitCollection* val)
{
    if (mMtdHits!= val) delete mMtdHits;
    mMtdHits = val;
}
//______________________________________________________________________________
void StMcEvent::setPxlHitCollection(StMcPxlHitCollection* val)
{
    if (mPxlHits!= val) delete mPxlHits;
    mPxlHits = val;
}   

//______________________________________________________________________________
void StMcEvent::setIstHitCollection(StMcIstHitCollection* val)
{
    if (mIstHits!= val) delete mIstHits;
    mIstHits = val;
}
   

//______________________________________________________________________________
void StMcEvent::setFgtHitCollection(StMcFgtHitCollection* val)
{
    if (mFgtHits!= val) delete mFgtHits;
    mFgtHits = val;
}
//______________________________________________________________________________
//______________________________________________________________________________
#define PrintHeader(Name,name) \
  const StMc ## Name ## HitCollection *name ## Coll = name ## HitCollection();\
  cout << "---------------------------------------------------------" << endl;\
  if (! all) {\
    cout << "StMc"#Name"HitCollection at " << (void*) name ## Coll             << endl; \
    cout << "Dumping collection size";                                  \
    cout << " and one hit only.";                                       \
    cout << endl;                                                       \
  }                                                                     \
  cout << "---------------------------------------------------------" << endl;\
  nhits = name ## Coll->numberOfHits();\
  cout << "# of hits in StMc"#Name"HitCollection = " << nhits << endl;
#define PrintHitCollection(Name,name)                                   \
  PrintHeader(Name,name)                                                \
    if ( name ## Coll &&  nhits) {                                      \
      nh = name ## Coll->hits().size();                                 \
      if (nh) {                                                         \
        if (! all) {nh = 1;  gotOneHit = kTRUE;}                        \
        for (i = 0; i < nh; i++) {                                      \
          if (! all) cout << *(name ## Coll->hits()[i]);                \
          else (name ## Coll->hits()[i])->Print();                      \
          cout << endl;                                                 \
          if (! all) {                                                  \
            cout << "Parent Track of this hit" << endl;                 \
            cout << *(name ## Coll->hits()[i]->parentTrack()) << endl;  \
          }                                                             \
        }                                                               \
      }                                                                 \
    }
#define PrintHitCollectionL(Name,name,layer,Layers)                     \
  PrintHeader(Name,name)                                                \
    if ( name ## Coll &&  nhits) {                                      \
      gotOneHit = kFALSE;                                               \
      for (k=0; !gotOneHit && k<name ## Coll->numberOf ## Layers(); k++) \
        if (name ## Coll->layer(k))                             \
          {                                                             \
            { nh = name ## Coll->layer(k)->hits().size();       \
              if (nh) {                                                 \
                  if (! all) {nh = 1;  gotOneHit = kTRUE;}              \
                  for (i = 0; i < nh; i++) {                            \
                    if (! all) cout << *(name ## Coll->layer(k)->hits()[i]); \
                    else   (name ## Coll->layer(k)->hits()[i])->Print(); \
                    cout << endl;                                       \
                    if (! all) {                                        \
                      cout << "Parent Track of this hit" << endl;       \
                      cout << *(name ## Coll->layer(k)->hits()[i]->parentTrack()) << endl; \
                    }                                                   \
                  }                                                     \
              }                                                         \
            }                                                           \
          }                                                             \
    }
#define PrintHitCollectionW(Name,name,ladder,Ladders,wafer,Wafers)      \
  PrintHeader(Name,name);                                               \
  if ( name ## Coll &&  nhits) {                                        \
    gotOneHit = kFALSE;                                                 \
    for (k=0; !gotOneHit && k<name ## Coll->numberOf ## Ladders(); k++) \
      if (name ## Coll->ladder(k))                                      \
        {                                                               \
          for (l = 0; gotOneHit && k<name ## Coll->ladder(k)->numberOf ## Wafers(); l++) \
            { nh = name ## Coll->ladder(k)->wafer(l)->hits().size();    \
              if (nh) {                                                 \
                if (! all) {nh = 1;  gotOneHit = kTRUE;}                \
                for (i = 0; i < nh; i++) {                              \
                  if (! all) cout << *(name ## Coll->ladder(k)->wafer(l)->hits()[i]); \
                  else   (name ## Coll->ladder(k)->wafer(l)->hits()[i])->Print(); \
                  cout << endl;                                         \
                  if (! all) {                                          \
                    cout << "Parent Track of this hit" << endl;         \
                    cout << *(name ## Coll->ladder(k)->wafer(l)->hits()[i]->parentTrack()) << endl; \
                  }                                                     \
                }                                                       \
              }                                                         \
            }                                                           \
        }                                                               \
  }

//________________________________________________________________________________
void StMcEvent::Print(Option_t *option) const {
  TString Opt(option);
  Int_t all = 0;
  if (Opt.Contains("all",TString::kIgnoreCase)) all = 1;
  //
  // Printing all the information of the components of StMcEvent.
  // For printing the simple data members of StMcEvent, use operator<<
  //
  //
  
  cout << "---------------------------------------------------------" << endl;
  cout << "StSPtrVecMcTrack"                                          << endl;
  if (! all) 
      cout << "Dumping first element in collection only (if available). " << endl;
  cout << "---------------------------------------------------------" << endl;
  Int_t Ntracks = tracks().size();
  cout << "collection size = " << Ntracks    << endl;
  if (Ntracks) {
      if (! all) { // Option for printing only the first track.
          cout << "---------------------------------------------------------" << endl;
          cout << "StMcTrack at " << tracks()[0] << endl;
          cout << "---------------------------------------------------------" << endl;
          cout << *(tracks()[0])                             << endl;
      }
      else { // Option for printing "all" tracks.
          tracks()[0]->Print("desc");
          for (int i = 0; i < Ntracks; i++) {
            (tracks()[i])->particleDefinition();
            (tracks()[i])->Print("all");
          }
      }
  }
  cout << "---------------------------------------------------------" << endl;
  if (! all) {
    cout << "StMcVertex" << endl;
    cout << "Dumping vertex info and first daughter track.            " << endl;
    cout << "---------------------------------------------------------" << endl;
    if (primaryVertex()) {
      cout << "Primary Vertex at " << primaryVertex() << endl;
      cout << "---------------------------------------------------------" << endl;
      cout << *primaryVertex() << endl;
      cout << "---------------------------------------------------------" << endl;
      cout << "First Daughter of Primary Vertex" << endl;
      cout << *(primaryVertex()->daughter(0)) << endl;
      cout << "---------------------------------------------------------" << endl;
    }
    else cout << "No Primary Vertex " << endl;
  }
  Int_t i1 = 0;
  Int_t nVertices = vertices().size();
  cout << "StSPtrVecMcVertex"                                         << endl;
  cout << "# of Vertices    : " << nVertices << endl;
  cout << "---------------------------------------------------------" << endl;
  if (nVertices > 1) {
    if (! all) {
      cout << "Daughters of second Vertex : " << vertices()[1]->numberOfDaughters() << endl;
      if (vertices()[1]->numberOfDaughters()) {
        cout << "First Daughter of this Vertex" << endl;
        cout << *(vertices()[1]->daughter(0)) << endl;
      }
      i1 = 1; 
      nVertices = 2;
      cout << "---------------------------------------------------------" << endl;
      cout << "Dumping vertices" << endl;
    }
    for (int i = i1; i < nVertices; i++) {
      if (! all) {
        cout << "---------------------------------------------------------" << endl;
        cout << "StMcVertex " << i+1 << " at " << vertices()[i]  << endl;
        cout << "---------------------------------------------------------" << endl;
        cout << *(vertices()[i]) << endl;
        cout << "---------------------------------------------------------" << endl;      
      }
      else {
        cout << "StMcVertex " << i+1 << " at ";
        (vertices()[i])->Print();
      }
    }
  }
  UInt_t       i, j, k, l, ii, nhits, nh;
  Bool_t             gotOneHit;
  PrintHeader(Tpc,tpc)
  if (tpcColl && nhits) {
      gotOneHit = kFALSE;
      for (k=0; !gotOneHit && k<tpcColl->numberOfSectors(); k++) {
          for (j=0; !gotOneHit && j<tpcColl->sector(k)->numberOfPadrows(); j++) {
              const StSPtrVecMcTpcHit &hits = tpcColl->sector(k)->padrow(j)->hits();
              nh = hits.size();
              if (nh) {
                  if (! all ) {
                      cout << *hits[0] << endl;
                      cout << "Parent Track of this hit" << endl;
                      cout << *(hits[0]->parentTrack()) << endl;
                      gotOneHit = kTRUE;
                      cout << "Dumping all the z coordinates and track key in this padrow" << endl;
                      cout << "Should be sorted according to z: " << endl;
                      cout << "---------------------------------------------------------" << endl;
                      for (i = 0; i < nh; i++) 
                          cout << "\t" << hits[i]->position().z() << ":" <<  hits[i]->parentTrack()->key();
                      cout << endl;
                  }
                  else for (i = 0; i < nh; i++) {
                    hits[i]->Print(); cout << endl;
                  }
              } // check for nh
              
          } // padrow loop
      }// sector loop
  }
  
  PrintHitCollectionL(Ftpc,ftpc,plane,Planes);
  PrintHitCollection(Rich,rich);
  
  PrintHeader(Svt,svt);
  nhits = svtColl->numberOfHits();
  if (svtColl && nhits) {
      
    gotOneHit = kFALSE;
    for (k=0; !gotOneHit && k<svtColl->numberOfBarrels(); k++)
        for (j=0; !gotOneHit && j<svtColl->barrel(k)->numberOfLadders(); j++)
            for (i=0; !gotOneHit && i<svtColl->barrel(k)->ladder(j)->numberOfWafers(); i++) {
                nh = svtColl->barrel(k)->ladder(j)->wafer(i)->hits().size();
                
                if (nh) {
                    if (! all) {nh = 1;  gotOneHit = kTRUE;}
                    for (ii = 0; ii < nh; ii++) {
                        if (! all) {
                          cout << *(svtColl->barrel(k)->ladder(j)->wafer(i)->hits()[ii]) << endl;
                          cout << "Parent track of this Hit" << endl;
                          cout << *(svtColl->barrel(k)->ladder(j)->wafer(i)->hits()[ii]->parentTrack()) << endl;
                        } else {(svtColl->barrel(k)->ladder(j)->wafer(i)->hits()[ii])->Print(); cout << endl;}
                    }
                }
            }
  }
  PrintHitCollectionW(Ssd,ssd,ladder,Ladders,wafer,Wafers);
  PrintHitCollection(Tof,tof);
  PrintHitCollection(BTof,btof);
  PrintHitCollection(Mtd,mtd);
  //PrintHitCollectionL(Pxl,pxl,sector,Sectors);
  PrintHitCollectionL(Ist,ist,layer,Layers);
  PrintHitCollectionL(Fgt,fgt,layer,Layers);
  PrintHitCollection(Etr,etr);
  
  TDataSet *mcEvent = (TDataSet *) this;
  TDataSetIter next(mcEvent);
  TDataSet *ds = 0;
  while( (ds = next())) {
    StMcEmcHitCollection *emcColl = dynamic_cast<StMcEmcHitCollection *>(ds);
    if (emcColl) emcColl->Print(option);
  }
}  
//______________________________________________________________________________
void StMcEvent::Browse(TBrowser *b)
{
  // Browse this event (called by TBrowser).
   for (auto jk=mBegColl; jk<mEndColl;jk++) { 
     TObject *obj = *jk; if (!obj) continue;
     b->Add(obj,obj->GetName());
   }
   TDataSet::Browse(b);
}
//______________________________________________________________________________
#undef  PrintHitCollection
#undef  PrintHitCollectionL
#undef  PrintHeader