STARLIGHT/starlight/TStarLight/TStarLight.cxx

   1 // -*- C++ -*-
   2 ////////////////////////////////////////////////////////////////////////
   3 //
   4 // Copyright 2013
   5 //
   6 ////////////////////////////////////////////////////////////////////////
   7 //
   8 // File and Version Information:
   9 // $Rev:: -1                      $: revision of last commit
  10 // $Authro::                      $: Author of last commit
  11 // $Date::                        $: Date of last commit
  12 //
  13 // Description:
  14 //     TStarLight.h is the include file defining data members and
  15 // functions specifically needed to implement an interface of STARlight
  16 // to ROOT's TGenerator (a part of ROOT's Virtual Monte Carlo.
  17 ////////////////////////////////////////////////////////////////////////
  18
  19 #include <string>
  20 #include <vector>
  21
  22 #include <Riostream.h>
  23 #include <TSystem.h>
  24 #include <TObjArray.h>
  25 #include <TClonesArray.h>
  26 #include <TParticle.h>
  27
  28 #include "inputParameters.h"
  29 #include "starlight.h"
  30 #include "TStarLight.h"
  31 #include "vector3.h"
  32
  33 ClassImp(TStarLight);
  34
  35 //----------------------------------------------------------------------
  36 TStarLight::TStarLight()
  37   : TGenerator()           // Default initlization of base class
  38   , fErrorStatus(0)        // Error status flag 0=OK
  39   , fConfigFileName("")    // Confiuration file name
  40   , fStarLight(NULL)       // STARlight simulation class.
  41   , fInputParameters(NULL) // Input to simulation class
  42   , fEvent() {}            // object holdng STARlight simulated event
  43
  44 //----------------------------------------------------------------------
  45 TStarLight::TStarLight(const char* name,         // The name of this object in the root name tables
  46                        const char* title,        // A title for this object
  47                        const char* slConfigFile) // file used to configure STARlight.
  48   : TGenerator(name, title)       // Default initlization of base class
  49   , fErrorStatus(0)               // Error status flag 0=OK
  50   , fConfigFileName(slConfigFile) // Confiuration file name
  51   , fStarLight(0)                 // STARlight simulation class.
  52   , fInputParameters(&inputParametersInstance) // Input to simulation class
  53   , fEvent()                      // object holding STARlight simulated event.
  54 {
  55   if (NULL == fInputParameters) {
  56     fErrorStatus = -5; // Init failed. Creating inputParamtere class
  57     Error("TStarLight", "creating inputParameters class failed");
  58     return;
  59   } // end if
  60   if (fConfigFileName != "" && gSystem->ExpandPathName(fConfigFileName)) { // if error
  61     fErrorStatus = -4; // Init failed. error in path name
  62     Error("TStarLight", "Error expanding path name='%s'", fConfigFileName.Data());
  63     return;
  64   } // end if
  65   fStarLight = new starlight();
  66   if (NULL == fStarLight) {
  67     fErrorStatus = -2; // Init failed. no simulation
  68     Error("TStarLight", "Failed to create simulator STARlight");
  69     return;
  70   } // end if
  71   if (fConfigFileName != "")
  72     fInputParameters->configureFromFile(slConfigFile);
  73 }
  74 //----------------------------------------------------------------------
  75 TStarLight::~TStarLight()
  76 {
  77   if (NULL != fStarLight) delete fStarLight;
  78   fStarLight = NULL;
  79   //   if (fInputParameters!=0) delete fInputParameters;
  80   //   fInputParameters = 0;
  81 }
  82
  83 //----------------------------------------------------------------------
  84 void TStarLight::GenerateEvent() {
  85
  86   if (NULL == fStarLight) {
  87     fErrorStatus = -1; // generate failed. No generator.
  88     Fatal("GenerateEvent", "TStarLight class/object not properly constructed");
  89     return;
  90   }
  91   fEvent = fStarLight->produceEvent();
  92 }
  93
  94 //----------------------------------------------------------------------
  95 Int_t TStarLight::ImportParticles(TClonesArray *part, // Pointer to array of particles to be filled.
  96                                   Option_t *opt) {    // A character array of options.
  97   // Return:
  98   //   The number of particles added to the TClonesArray *part.
  99
 100   if (NULL == part)
 101     return 0;
 102
 103   TClonesArray &clonesParticles = *part;
 104   clonesParticles.Clear();
 105
 106   Int_t nVtx(0);
 107   Double_t vtx(0), vty(0), vtz(0), vtt(0);
 108   const std::vector<vector3>* slVtx = fEvent.getVertices();
 109   if (NULL == slVtx) { // not vertex assume 0,0,0,0;
 110     vtx = vty = vtz = vtt = 0.0;
 111   } else { // a vertex exits
 112     slVtx = fEvent.getVertices();
 113     nVtx  = slVtx->size();
 114   }
 115   const std::vector<starlightParticle>* slPartArr = fEvent.getParticles();
 116   const Int_t npart(slPartArr->size());
 117   if (!strcmp(opt,"") || !strcmp(opt,"Final")) {
 118     for (Int_t ipart=0; ipart<npart; ++ipart) {
 119       const starlightParticle* slPart = &(slPartArr->at(ipart));
 120       if (Stable(slPart->getPdgCode())) {
 121         //slPart->setParent(-1);
 122         if (nVtx<1) { // No verticies
 123           vtx = vty = vtz = vtt = 0.0;
 124         } else {
 125           vtx = (slVtx->at(ipart <nVtx ? ipart : 0)).X();
 126           vty = (slVtx->at(ipart <nVtx ? ipart : 0)).Y();
 127           vtz = (slVtx->at(ipart <nVtx ? ipart : 0)).Z();
 128           vtt = 0.0; // no time given.
 129         } // end if
 130         new(clonesParticles[ipart]) TParticle(slPart->getPdgCode(),
 131                                               slPart->getStatus(),
 132                                               -1,/*slPart->getParent(),*/
 133                                               -1,
 134                                               slPart->getFirstDaughter(),
 135                                               slPart->getLastDaughter(),
 136                                               slPart->GetPx(),
 137                                               slPart->GetPy(),
 138                                               slPart->GetPz(),
 139                                               slPart->GetE(),
 140                                               vtx,vty,vtz,vtt);
 141       }
 142     }
 143   } else if (!strcmp(opt,"ALL")) { // store all particles
 144     for (Int_t ipart=0; ipart<npart; ++ipart) {
 145       const starlightParticle* slPart = &(slPartArr->at(ipart));
 146       //slPart->setParent(-1);
 147       if (nVtx < 1) { // No verticies
 148         vtx = vty = vtz = vtt = 0.0;
 149       } else {
 150         vtx = (slVtx->at((ipart<nVtx?ipart:0))).X();
 151         vty = (slVtx->at((ipart<nVtx?ipart:0))).Y();
 152         vtz = (slVtx->at((ipart<nVtx?ipart:0))).Z();
 153         vtt = 0.0; // no time given.
 154       }
 155       new(clonesParticles[ipart]) TParticle(slPart->getPdgCode(),
 156                                             slPart->getStatus(),
 157                                             -1,/*slPart->getParent(),*/
 158                                             -1,
 159                                             slPart->getFirstDaughter(),
 160                                             slPart->getLastDaughter(),
 161                                             slPart->GetPx(),
 162                                             slPart->GetPy(),
 163                                             slPart->GetPz(),
 164                                             slPart->GetE(),
 165                                             vtx,vty,vtz,vtt);
 166     } // end for ipart
 167   } else {
 168     fErrorStatus = -1; // Import particles failed unknown option
 169     Error("ImportParticles", "Unknown option '%s'", opt);
 170   } // end if opt
 171   return npart;
 172 }
 173
 174 //----------------------------------------------------------------------
 175 TObjArray* TStarLight::ImportParticles(Option_t *opt) { // A character array of options
 176   Int_t nVtx(0);
 177   Double_t vtx(0), vty(0), vtz(0), vtt(0);
 178   const std::vector<vector3>* slVtx(fEvent.getVertices());
 179   if (slVtx == 0) { // not vertex assume 0,0,0,0;
 180     vtx = vty = vtz = vtt = 0.0;
 181   } else { // a vertex exits
 182     slVtx = fEvent.getVertices();
 183     nVtx = slVtx->size();
 184   } // end if
 185   const std::vector<starlightParticle>* slPartArr(fEvent.getParticles());
 186   const Int_t npart(fEvent.getParticles()->size());
 187   if (!strcmp(opt, "") || !strcmp(opt, "Final")) {
 188     for (Int_t ipart=0; ipart<npart; ipart++) {
 189       const starlightParticle* slPart(&(slPartArr->at(ipart)));
 190       if (Stable(slPart->getPdgCode())) {
 191         //slPart->setParent(-1);
 192         if (nVtx < 1) { // No verticies
 193           vtx = vty = vtz = vtt = 0.0;
 194         } else {
 195           vtx = (slVtx->at((ipart < nVtx ? ipart : 0))).X();
 196           vty = (slVtx->at((ipart < nVtx ? ipart : 0))).Y();
 197           vtz = (slVtx->at((ipart < nVtx ? ipart : 0))).Z();
 198           vtt = 0.0; // no time given.
 199         } // end if
 200         TParticle *p = new TParticle(slPart->getPdgCode(),
 201                                      slPart->getStatus(),
 202                                      -1,/*slPart->getParent(),*/
 203                                      -1,
 204                                      slPart->getFirstDaughter(),
 205                                      slPart->getLastDaughter(),
 206                                      slPart->GetPx(),
 207                                      slPart->GetPy(),
 208                                      slPart->GetPz(),
 209                                      slPart->GetE(),
 210                                      vtx,vty,vtz,vtt);
 211         fParticles->Add(p);
 212       }
 213     } // end for ipart
 214   } else if (!strcmp(opt,"ALL")) {// store all particles
 215     for(Int_t ipart=0;ipart<npart;ipart++) {
 216       const starlightParticle* slPart(&(slPartArr->at(ipart)));
 217       //slPart->setParent(-1);
 218       if (nVtx < 1) { // No verticies
 219         vtx = vty = vtz = vtt = 0.0;
 220       } else {
 221         vtx = (slVtx->at((ipart < nVtx ? ipart : 0))).X();
 222         vty = (slVtx->at((ipart < nVtx ? ipart : 0))).Y();
 223         vtz = (slVtx->at((ipart < nVtx ? ipart : 0))).Z();
 224         vtt = 0.0; // no time given.
 225       } // end if
 226       TParticle* p = new TParticle(slPart->getPdgCode(),
 227                                    slPart->getStatus(),
 228                                    -1,/*slPart->getParent(),*/
 229                                    -1,
 230                                    slPart->getFirstDaughter(),
 231                                    slPart->getLastDaughter(),
 232                                    slPart->GetPx(),
 233                                    slPart->GetPy(),
 234                                    slPart->GetPz(),
 235                                    slPart->GetE(),
 236                                    vtx,vty,vtz,vtt);
 237       fParticles->Add(p);
 238     } // end for ipart
 239   }else{
 240     fErrorStatus = -1; // Import particles failed unknown option
 241     Error("ImportParticles", "Unknown option '%s'", opt);
 242   } // end if opt
 243   return fParticles;
 244 }
 245 //----------------------------------------------------------------------
 246 void TStarLight::SetParameter(const char* line) {
 247   const std::string sl(line);
 248   if (not fInputParameters->setParameter(sl))
 249     Fatal("SetParameter", "cannot set parameter: '%s'", line);
 250 }
 251 void TStarLight::SetParameter(const char* name,
 252                               Double_t val) {
 253   const TString line(TString::Format("%s = %e", name, val));
 254   const std::string sl(line.Data());
 255   if (not fInputParameters->setParameter(sl))
 256     Fatal("SetParameter", "cannot set parameter: '%s'", name);
 257 }
 258
 259 //----------------------------------------------------------------------
 260 Double_t TStarLight::GetParameter(const char* name) const {
 261   const std::string sn(name);
 262   std::stringstream ioss;
 263   fInputParameters->write(ioss);
 264   std::string line;
 265   while (getline(ioss, line)) {
 266     const size_t index_colon(line.find(':'));
 267     if (index_colon == std::string::npos)
 268       continue;
 269     std::string key(line.substr(0,index_colon));
 270     if (key != sn) continue;
 271     std::string val(line.substr(index_colon+1, std::string::npos));
 272     std::istringstream iss(val);
 273     Double_t result(0);
 274     iss >> result;
 275     return result;
 276   }
 277   Fatal("GetParameter", "parameter '%s' not found", name);
 278   return 0.0;
 279 }