loon: RerootToTruthModule.cxx Source File

00001 
00002 // $Id: RerootToTruthModule.cxx,v 1.51 2008/05/16 15:50:29 rhatcher Exp $
00003 //
00004 // A JobControl Module for filling Truth from REROOT
00005 //
00006 // rhatcher@fnal.gov
00008 
00009 #include <cassert>
00010 
00011 #include <RerootExodus/RerootToTruthModule.h>
00012 #include <RerootExodus/RerootExodus.h>
00013 
00014 #include <MinosObjectMap/MomNavigator.h>
00015 
00016 #include <MessageService/MsgService.h>
00017 #include <JobControl/JobCModuleRegistry.h>
00018 
00019 #include <Conventions/Mphysical.h>
00020 #ifdef EVENT_KINEMATICS_PKGS
00021 #include <EventKinematics/NuEvtKin.h>
00022 #include <EventKinematics/StdHepUtil.h>
00023 using namespace StdHepUtil;
00024 #endif
00025 #include <REROOT_Classes/REROOT_StdHep.h>
00026 #include <REROOT_Classes/REROOT_StdHepHead.h>
00027 #include <REROOT_Classes/REROOT_NeuKin.h>
00028 #include <REROOT_Classes/REROOT_FluxInfo.h>
00029 #include <REROOT_Classes/REROOT_FluxWgt.h>
00030 #include <Record/SimSnarlRecord.h>
00031 #include <Record/SimSnarlHeader.h>
00032 #include <UgliGeometry/UgliGeomHandle.h>
00033 #include <UgliGeometry/UgliStripHandle.h>
00034 
00035 #include <Digitization/DigiScintHit.h>
00036 #include <Digitization/DigiRerootInfo.h>
00037 
00038 #include <DatabaseInterface/DbiTableProxyRegistry.h>
00039 
00040 #include <Util/LoadMinosPDG.h>
00041 #include <Util/UtilPDG.h>
00042 #include <Util/UtilHepevt.h>
00043 
00044 #include <TMath.h>
00045 #include <TString.h>
00046 #include <TSystem.h>
00047 #include <TDatabasePDG.h>
00048 #include <TParticle.h>
00049 #include <TLorentzVector.h>
00050 #include <TROOT.h>
00051 #include <TApplication.h>
00052 
00053 ClassImp(RerootToTruthModule)
00054 
00055 //......................................................................
00056 
00057 CVSID("$Id: RerootToTruthModule.cxx,v 1.51 2008/05/16 15:50:29 rhatcher Exp $");
00058 JOBMODULE(RerootToTruthModule, "RerootToTruthModule",
00059          "Builds RawRecord with RawDigitDataBlock");
00060 
00061 // by default no override function
00062 static std::string fgGlobalNewVtxFunction = "";
00063 std::string RerootToTruthModule::GetGlobalNewVtxFunction()
00064 { return fgGlobalNewVtxFunction; }
00065 void RerootToTruthModule::SetGlobalNewVtxFunction(std::string funcname)
00066 { fgGlobalNewVtxFunction = funcname; }
00067 
00068 // default for fMakeDigiScintList
00069 static int fgGlobalMakeDigiScintList = 1;
00070 int RerootToTruthModule::GetGlobalMakeDigiScintList()
00071 { return fgGlobalMakeDigiScintList; }
00072 void RerootToTruthModule::SetGlobalMakeDigiScintList(int flag)
00073 { fgGlobalMakeDigiScintList = flag; }
00074 
00075 
00076 //......................................................................
00077 
00078 RerootToTruthModule::RerootToTruthModule() 
00079     : detector(Detector::kUnknown)
00080     , fGetApplyRollback(true)
00081     , fAnaPrintStdHep(true)
00082     , fAnaPrintStdHepHead(false)
00083     , fAnaPrintNeuKin(true)
00084     , fAnaPrintNuEvtKin(true)
00085     , fAnaPrintFluxInfo(true)
00086     , fAnaPrintFluxWgt(true)
00087     , fAnaPrintScintHit(false)
00088     , fIonScheme(UtilPDG::kIonUnchanged)
00089     , fFixCharmDKStatus(true)
00090     , fDropGeantEntries(false)
00091     , fDropTauDecayProducts(false)
00092     , fDropCharmDecayProducts(false)
00093     , fDropStatus999(false)
00094     , fStdHepEditString("")
00095     , fMakeDigiScintHitList(fgGlobalMakeDigiScintList)
00096     , fNewVtxFunction(fgGlobalNewVtxFunction)
00097 {
00098    // construct a new "RerootToTruthModule" JobControl module
00099 
00100    // initialize detector to non-legal value
00101 
00102    LoadMinosPDG();
00103 
00104 }
00105 
00106 //......................................................................
00107 
00108 RerootToTruthModule::~RerootToTruthModule() 
00109 {
00110   MSG("Exodus", Msg::kVerbose) << "RerootToTruthModule::Destructor\n";
00111 
00112 }
00113 
00114 //......................................................................
00115 // Helper function, creates a TClonesArray of TParticles from StdHep
00116 // list. 
00117 static TClonesArray* make_std_hep_list() 
00118 {
00119    // Produce the equivalent of StdHep
00120    const TClonesArray *stdheps = RerootExodus::GetStdHepList();
00121    REROOT_StdHep *stdhep;
00122    Int_t nstdheps = stdheps->GetLast()+1;
00123 
00124    // New TClonesArray of the right size
00125    TClonesArray* newstdhep = new TClonesArray("TParticle",nstdheps);
00126    newstdhep->SetName("StdHep");
00127    Int_t entry = 0;
00128 
00129    for (Int_t ihep=0; ihep<nstdheps; ihep++) {
00130       stdhep = dynamic_cast<REROOT_StdHep*>(stdheps->UncheckedAt(ihep));
00131 
00132       const Int_t ipdg      = UtilPDG::stdIonNumber(stdhep->IdHEP());
00133       const Int_t status    = stdhep->IstHEP();
00134       const Int_t mother1   = stdhep->FirstParent();
00135       const Int_t mother2   = stdhep->LastParent();
00136       const Int_t daughter1 = stdhep->FirstChild();
00137       const Int_t daughter2 = stdhep->LastChild();
00138       const Double_t px     = stdhep->Px()   *Munits::GeV;
00139       const Double_t py     = stdhep->Py()   *Munits::GeV;
00140       const Double_t pz     = stdhep->Pz()   *Munits::GeV;
00141     //const Double_t mass   = stdhep->Mass() *Munits::GeV;
00142       const Double_t etot   = stdhep->E()    *Munits::GeV; 
00143     //const Double_t ecalc  = TMath::Sqrt(px*px+py*py+pz*pz+mass*mass);
00144       const Double_t vx     = stdhep->Xmm()  *Munits::mm;
00145       const Double_t vy     = stdhep->Ymm()  *Munits::mm;
00146       const Double_t vz     = stdhep->Zmm()  *Munits::mm;
00147       // the time in StdHep was ct in "mm", convert to base time units
00148       const Double_t time   = stdhep->Tprod() *Munits::mm / Mphysical::c_light;
00149 
00150       // new with placement
00151       new((*newstdhep)[entry++])
00152          TParticle(ipdg,status,mother1,mother2,daughter1,daughter2,
00153                    px,py,pz,etot,vx,vy,vz,time);
00154 
00155    }
00156    return newstdhep;
00157 }
00158 
00159 static TClonesArray* make_neukin_list() 
00160 {
00161    // Produce the TClonesArray of REROOT_NeuKin objects
00162    const TClonesArray *neukins = RerootExodus::GetNeuKinList();
00163    Int_t nneukins = neukins->GetLast()+1;
00164 
00165    // New TClonesArray of the right size
00166    TClonesArray* newneukins = new TClonesArray("REROOT_NeuKin",nneukins);
00167    newneukins->SetName("NeuKinList");
00168    Int_t entry = 0;
00169 
00170    for (Int_t ikin=0; ikin<nneukins; ikin++) {
00171       const REROOT_NeuKin* old_neukin = 
00172         dynamic_cast<const REROOT_NeuKin*>(neukins->UncheckedAt(ikin));
00173       // new with placement
00174       new((*newneukins)[entry++])REROOT_NeuKin(*old_neukin);
00175    }
00176    return newneukins;
00177 }
00178 
00179 static TClonesArray* make_fluxinfo_list() 
00180 {
00181    // Produce the TClonesArray of REROOT_FluxInfo objects
00182    const TClonesArray *fluxinfos = RerootExodus::GetFluxInfoList();
00183    Int_t nfluxinfos = fluxinfos->GetLast()+1;
00184 
00185    // New TClonesArray of the right size
00186    TClonesArray* newfluxinfos = new TClonesArray("REROOT_FluxInfo",nfluxinfos);
00187    newfluxinfos->SetName("FluxInfoList");
00188    Int_t entry = 0;
00189 
00190    for (Int_t ifi=0; ifi<nfluxinfos; ifi++) {
00191       const REROOT_FluxInfo* old_fluxinfo = 
00192         dynamic_cast<const REROOT_FluxInfo*>(fluxinfos->UncheckedAt(ifi));
00193       // new with placement
00194       new((*newfluxinfos)[entry++])REROOT_FluxInfo(*old_fluxinfo);
00195    }
00196    return newfluxinfos;
00197 }
00198 
00199 static TClonesArray* make_fluxwgt_list() 
00200 {
00201    // Produce the TClonesArray of REROOT_FluxWgt objects
00202    const TClonesArray *fluxwgts = RerootExodus::GetFluxWgtList();
00203    Int_t nfluxwgts = fluxwgts->GetLast()+1;
00204 
00205    // New TClonesArray of the right size
00206    TClonesArray* newfluxwgts = new TClonesArray("REROOT_FluxWgt",nfluxwgts);
00207    newfluxwgts->SetName("FluxWgtList");
00208    Int_t entry = 0;
00209 
00210    for (Int_t ifw=0; ifw<nfluxwgts; ifw++) {
00211       const REROOT_FluxWgt* old_fluxwgt = 
00212         dynamic_cast<const REROOT_FluxWgt*>(fluxwgts->UncheckedAt(ifw));
00213       // new with placement
00214       new((*newfluxwgts)[entry++])REROOT_FluxWgt(*old_fluxwgt);
00215    }
00216    return newfluxwgts;
00217 }
00218 
00219 #ifdef EVENT_KINEMATICS_PKGS
00220 Int_t next_doc_nu(TClonesArray* stdhep, Int_t ihep_current)
00221 {
00222   // find the next isDocStatus, isNeutrino entry in StdHep
00223   // *AFTER* the current location, or one beyond end of array
00224   Int_t last = stdhep->GetLast();
00225   for (int ih = ihep_current+1; ih <= last; ++ih) {
00226     const TParticle* part = (const TParticle*)stdhep->UncheckedAt(ih);
00227     if ( isDocStatus(part) && isNeutrino(part) ) return ih;
00228   }
00229   return last+1;
00230 }
00231 
00232 static TClonesArray* make_NuEvtKin_list(TClonesArray* stdhep) 
00233 {
00234    // Produce the TClonesArray of NuEvtKin objects
00235   
00236    // verify that StdHep is "StdHep"
00237    TString stdhep_name = stdhep->GetName();
00238    if (stdhep_name != "StdHep") {
00239      MSG("Exodus",Msg::kWarning)
00240        << "make_NuEvtKin_list was passed TClonesArray named \""
00241        << stdhep_name << "\" rather than \"StdHep\"" << endl;
00242      return 0;
00243    }
00244      
00245    const TClonesArray *neukins = RerootExodus::GetNeuKinList();
00246    Int_t nneukins = neukins->GetLast()+1;
00247 
00248    // New TClonesArray of the right size
00249    TClonesArray* nuevtkins = new TClonesArray("NuEvtKin",nneukins);
00250    nuevtkins->SetName("NuEvtKinList");
00251    Int_t entry = 0;
00252    Int_t ihep  = 0;
00253 
00254    for (Int_t ikin=0; ikin<nneukins; ikin++) {
00255       const REROOT_NeuKin* old_neukin = 
00256         dynamic_cast<const REROOT_NeuKin*>(neukins->UncheckedAt(ikin));
00257       if ( ! old_neukin ) continue;
00258 
00259       NuEvtKin::Scatter_t scatter = NuEvtKin::kUnknownScatter;
00260       switch (old_neukin->IResonance()) {  // translate Neugen3 codes
00261       case 1001: scatter = NuEvtKin::kQEScatter;         break;
00262       case 1002: scatter = NuEvtKin::kResonantPiScatter; break;
00263       case 1003: scatter = NuEvtKin::kDIScatter;         break;
00264       case 1004: scatter = NuEvtKin::kCoherentPiScatter; break;
00265       }
00266 
00267       NuEvtKin::Interaction_t inter = NuEvtKin::kUnknownInteraction;
00268       switch (old_neukin->IAction()) {
00269       case 0: inter = NuEvtKin::kWeakNC; break;
00270       case 1: inter = NuEvtKin::kWeakCC; break;
00271       }
00272 
00273       //Int_t ipdgNu      = old_neukin->INu();
00274       Int_t ipdgNuNoOsc = old_neukin->INuNoOsc();
00275       Float_t x         = old_neukin->X();
00276       Float_t y         = old_neukin->Y();
00277       Float_t Q2        = old_neukin->Q2();
00278       Float_t W2        = old_neukin->W2();
00279 
00280       // new with placement
00281       new((*nuevtkins)[entry])NuEvtKin(stdhep,ipdgNuNoOsc,
00282                                        scatter,inter,x,y,Q2,W2);
00283 
00284       NuEvtKin* currentNuEvtKin = (NuEvtKin*)((*nuevtkins)[entry]);
00285       // now fix up the entry by telling it what StdHep entries to use
00286       // assume segments of StdHep start with the documentation neutrino entry
00287       Int_t ihep_next_nu = next_doc_nu(stdhep,ihep);
00288       for ( ; ihep < ihep_next_nu ; ++ihep ) 
00289         currentNuEvtKin->LinkStdHepIndexToEvent(ihep);
00290 
00291       // prepare for next entry
00292       entry++;
00293    }
00294    return nuevtkins;
00295 }
00296 #endif
00297 
00298 static double flshit_path_length(REROOT_FLSHit* h)
00299 {                               // assuming a strait track!
00300     double dx = h->XEnd() - h->XBegin();
00301     double dy = h->YEnd() - h->YBegin();
00302     double dz = h->ZEnd() - h->ZBegin();
00303     return sqrt(dx*dx+dy*dy+dz*dz)*Munits::cm;
00304 }
00305 static double flshit_track_energy(REROOT_FLSHit* h)
00306 {
00307     double p = h->Ptot()*Munits::GeV;
00308     TDatabasePDG *pdg = TDatabasePDG::Instance();
00309     int ipdg = h->IPDG();
00310     int jjj = ipdg%1000;
00311 
00312     // silently fix bad value due to real<->integer conversion in
00313     // original FORTRAN MC 
00314     if (jjj and ipdg > 100000000) 
00315       ipdg += ((jjj<500) ? (-jjj) : (1000-jjj));
00316 
00317     TParticlePDG *part = pdg->GetParticle(ipdg);
00318     if (!part) {
00319         static int nmsg = 25;
00320         if (nmsg) {
00321           MSG("Exodus", Msg::kWarning) 
00322             << "unknown particle with PDG #" << h->IPDG() 
00323             << " (" << ipdg << ")"
00324             << endl;
00325           MSG("Exodus", Msg::kWarning) << "Assuming zero energy!!!\n";
00326           if (--nmsg == 0) 
00327             MSG("Exodus", Msg::kWarning)
00328               << "  ... last warning of this type" << endl;
00329         }
00330         return 0.0;
00331     }
00332     double m = part->Mass()*Munits::GeV;
00333     return sqrt(p*p+m*m);
00334 }
00335 static void flshit_massage_local(PlexStripEndId seid, float x, float y, float z,
00336                                  float &xout, float &yout, float &zout)
00337 {
00338     // A painful little dance....
00339 
00340     // near detector REROOT data has phantom strips
00341     // that will have no UgliStripHandle ... 
00342     if ( seid.GetStrip() >= seid.GetNumStrips() ) {
00343       static int nmsg = 10;
00344       if ( nmsg > 0 ) {
00345         MSG("Exodus",Msg::kWarning)
00346           << "flshit_massage_local called for phantom strip "
00347           << seid << endl;
00348         nmsg--;
00349         if (nmsg == 0) 
00350           MSG("Exodus",Msg::kWarning) 
00351             << " ... last warning of this type" << endl;
00352       }
00353       // can't do anything sensible so leave values unmodified
00354       xout = x; yout = y; zout = z;
00355       return;
00356     }
00357 
00358     // starting with local REROOT position, go to global Ugli
00359     TVector3 local (x,y,z);
00360     TVector3 global = RerootExodus::SEIdLocalToGlobal(seid,local);
00361 
00362     global[0] *= Munits::cm;
00363     global[1] *= Munits::cm;
00364     global[2] *= Munits::cm;
00365 
00366     // Then from global Ugli to local Ugli
00367     VldContext vld = RerootExodus::BuildVldContext();
00368     UgliGeomHandle ugh(vld);
00369     assert (ugh.IsValid());
00370     UgliStripHandle ush = ugh.GetStripHandle(seid);
00371     local = ush.GlobalToLocal(global);
00372 
00373     xout = local[0];
00374     yout = local[1];
00375     zout = local[2];
00376 }
00377 
00378 static TClonesArray* make_digi_scint_hit_list()
00379 {
00380     const TClonesArray *old_arr = RerootExodus::GetFLSHitList();
00381     int size = old_arr->GetLast()+1;
00382 
00383     TClonesArray *new_arr = new TClonesArray("DigiScintHit",size);
00384     new_arr->SetName("DigiScintHits");
00385 
00386     for (int ind = 0; ind < size; ++ind) {
00387         REROOT_FLSHit* flshit = 
00388             dynamic_cast<REROOT_FLSHit*>(old_arr->UncheckedAt(ind));
00389         if (!flshit) continue;
00390 
00391         // IAB=1 is generally West side except old 'F' plane 'V' view case
00392         PlexStripEndId seid = 
00393           RerootExodus::PECAB2SEId(flshit->IPln(),
00394                                    flshit->IExtr(),
00395                                    flshit->ICell(),1);  
00396 
00397         // force StripEnd::kWest; we don't really care about the
00398         // actual end but PlexStripEndId::BuildPlnStripEndKey() is
00399         // being used later and it complains about use of NearDet 
00400         // East ends (which don't exist for readout purposes).
00401         seid.SetEnd(StripEnd::kWest);
00402 
00403         double ds = flshit_path_length(flshit);
00404         double e = flshit_track_energy(flshit);
00405 
00406         float x1,y1,z1,x2,y2,z2;
00407 
00408         // determine the plane type
00409         char plnkey = RerootExodus::PlaneName(seid)[0];
00410 
00411         if ( plnkey == 'M' ) {
00412           // just need to convert from GMINOS units to standard units.
00413           x1 = flshit->XBegin() * Munits::cm;
00414           y1 = flshit->YBegin() * Munits::cm;
00415           z1 = flshit->ZBegin() * Munits::cm;
00416           x2 = flshit->XEnd() * Munits::cm;
00417           y2 = flshit->YEnd() * Munits::cm;
00418           z2 = flshit->ZEnd() * Munits::cm;
00419         }
00420         else {
00421           // RWH: 031219
00422           // I don't know why there is this massaging going on
00423           // but one certainly doesn't want it for GMINOS 'M' planes
00424           // where {X|Y|Z}{Begin|End} are actually local coords
00425           flshit_massage_local(seid,flshit->XBegin(),flshit->YBegin(),flshit->ZBegin(),
00426                                x1,y1,z1);
00427           flshit_massage_local(seid,flshit->XEnd(),flshit->YEnd(),flshit->ZEnd(),
00428                                x2,y2,z2);
00429         }
00430 
00431         new ((*new_arr)[ind]) DigiScintHit(flshit->ITrack(),
00432                                            flshit->IPDG(),
00433                                            e,
00434                                            seid,
00435                                            flshit->TOFG(),
00436                                            x1,y1,z1,
00437                                            flshit->TOFG() + ds/Munits::c_light,
00438                                            x2,y2,z2,
00439                                            ds,
00440                                            flshit->ELoss());
00441     }
00442 
00443     if (new_arr->GetLast() + 1 != size) {
00444         MSG("Exodus", Msg::kWarning) << 
00445             "make_digi_scint_hit_list: some objects not convertable to FLSHits\n";
00446     }
00447         
00448     return new_arr;
00449 }
00450 
00451 static void print_tclones_stdhep(const TClonesArray* newstdhep)
00452 {
00453    // Loop over new array ... print each element
00454    Bool_t doheader = kTRUE;
00455    int nhep = newstdhep->GetLast() + 1;
00456    printf(" nhep = %d\n",nhep);
00457    for (Int_t i=0; i<nhep; i++) {
00458       TParticle* apart = dynamic_cast<TParticle*>(newstdhep->At(i));
00459       if (apart) {
00460          //rwh: does a cruddy job: apart->Print();
00461 
00462          static Char_t alt[64];
00463 
00464          Int_t status = apart->GetStatusCode();
00465          Int_t ipdg   = apart->GetPdgCode();
00466          const Char_t* name = apart->GetName();
00467          if (name[0]=='X' && name[1]=='X' && name[2] == 'X') {
00468             sprintf(alt,"%d",ipdg);
00469             name = alt;
00470          }
00471          Int_t first_parent = apart->GetMother(0);
00472          Int_t last_parent  = apart->GetMother(1);
00473          // PTSim uses convention of 2nd parent as -1 if there is only one
00474          if (first_parent!=last_parent && last_parent!=-1) 
00475            printf(" next line had different parents %d %d:\n",
00476                   first_parent,last_parent);
00477          Int_t first_child  = apart->GetDaughter(0);
00478          Int_t last_child   = apart->GetDaughter(1);
00479          Double_t px = apart->Px();
00480          Double_t py = apart->Py();
00481          Double_t pz = apart->Pz();
00482          Double_t energy = apart->Energy();
00483          //Double_t calcmass = apart->GetCalcMass();
00484          //Double_t mass = (apart->GetPDG())?apart->GetMass():-99999.;
00485          Double_t vx = apart->Vx();
00486          Double_t vy = apart->Vy();
00487          Double_t vz = apart->Vz();
00488          Double_t tprod = apart->T();
00489 
00490          //if (doheader) {
00491          //   doheader = kFALSE;
00492          //   //      12345678901234567890123456789012345678901234567890123456789012345678901234567890
00493          //   //                           1234123412341234 123456789A 123456789A 123456789A 12345678
00494          //   printf("ihep stat type           [parents]         px          py          pz     energy\n");
00495          //   printf("          ipdg           [childrn]         vx          vy          vz      tprod\n");
00496          //}
00497          //printf("%4d(%3d) %-10.10s     %4d%4d",
00498          //       i,status,name,first_parent,last_parent);
00499          //printf(" %11.4g %11.4g %11.4g%11.4g\n",px,py,pz,energy);
00500          //printf("          %10d     %4d%4d",ipdg,first_child,last_child);
00501          //printf(" %11.4g %11.4g %11.4g%11.4g\n",vx,vy,vz,tprod);
00502 
00503          if (doheader) {
00504             doheader = kFALSE;
00505             //      12345678901234567890123456789012345678901234567890123456789012345678901234567890
00506             //                           1234123412341234 123456789A 123456789A 123456789A 12345678
00507             printf("ihep stat type        <parents  >          px          py          pz     energy\n");
00508             printf("                      [children ]          vx          vy          vz      tprod\n");
00509             //                 <1234 1234> [1234 1235]
00510          }
00511          printf("%4d(%4d)%11.11s <%4d %4d>",
00512                 i,status,name,first_parent,last_parent);
00513          printf(" %11.6g %11.6g %11.6g%11.6g\n",px,py,pz,energy);
00514          printf("          %11d [%4d %4d]",
00515                 ipdg,first_child,last_child);
00516          printf(" %11.6g %11.6g %11.6g %10.4e\n",vx,vy,vz,tprod);
00517 
00518 
00519       } else {
00520          MSG("Exodus", Msg::kInfo) << "RerootToTruthModule new StdHep element " << i << " is empty" << endl;
00521       }
00522    }
00523 }
00524 
00525 static void print_tclones_stdhephead() 
00526 {
00527    // Produce the TClonesArray of REROOT_FluxWgt objects
00528    const TClonesArray *stdhepheads = RerootExodus::GetStdHepHeadList();
00529    Int_t nstdhepheads = stdhepheads->GetLast()+1;
00530 
00531    for (Int_t ishh=0; ishh<nstdhepheads; ishh++) {
00532       const REROOT_StdHepHead* old_shh = 
00533         dynamic_cast<const REROOT_StdHepHead*>(stdhepheads->UncheckedAt(ishh));
00534       cout << " entry " << setw(3) << ishh 
00535            << " ID= " << setw(3) << old_shh->ID()
00536            << " NevHEP= " << setw(7) << old_shh->NevHEP()
00537            << " NHEP= " << setw(5) << old_shh->NHEP()
00538            << endl;
00539    }
00540 
00541 }
00542 //......................................................................
00543 
00544 JobCResult RerootToTruthModule::Get(MomNavigator *mom)
00545 {
00546    // Create SimSnarlRecord from the current REROOT event
00547    // available via gMINFast.
00548 
00549    detector = RerootExodus::GetDetector();
00550 
00551    // For a PlexHandle we need a context
00552    VldContext vldc = RerootExodus::BuildVldContext();
00553 
00554    Int_t run      = RerootExodus::GetRunNo();
00555    Short_t subrun = RerootExodus::GetSubRunNum();
00556    Int_t snarl    = RerootExodus::GetSnarlNum(); //tie to tree entry
00557 
00558    MsgService::Instance()->SetCurrentRunSnarl(run,snarl);
00559 
00560    // produce a SimSnarlRecord
00561    // have it adopt (secure) what we've produed
00562   
00563    Short_t runtype = 0;
00564    UInt_t  trigsrc = 0;
00565    UInt_t  errcode = 0;
00566    Int_t   timeframe = RerootExodus::GetTimeFrame(vldc);
00567    Int_t   spilltype = -1;
00568 
00570    // Find or make the SimSnarl
00571    
00572    SimSnarlRecord* record = 0;
00573    
00574    // See if a SimSnarl already exists. 
00575    TObject* tobj;
00576    TIter    fragiter = mom->FragmentIter();
00577    while( ( tobj = fragiter.Next() ) ) {
00578      record = dynamic_cast<SimSnarlRecord*>(tobj);
00579      if(record) break;
00580    }
00581    
00582   // If not, make one.
00583   if(record == 0) {
00584     VldTimeStamp mcGenTime(RerootExodus::GetLastEventHistoryTimeStamp());
00585     std::string  mcCodename(RerootExodus::GetGminosCodeName().Data());
00586     std::string  mcHostname(RerootExodus::GetGminosHostName().Data());
00587     SimSnarlHeader simheader(vldc,run,subrun,runtype,errcode,
00588                              snarl,trigsrc,timeframe,spilltype,
00589                              mcGenTime,mcCodename,mcHostname);
00590     record = new SimSnarlRecord(simheader);
00591     // Add the GMINOS information to the record
00592     RecJobHistory& jobhist 
00593                     = const_cast<RecJobHistory&>(record->GetJobHistory());
00594     jobhist.CreateJobRecord(RecJobHistory::kGMinos,mcCodename.c_str(),"",
00595                             mcHostname.c_str(),mcGenTime);
00596     
00597     MSG("Exodus",Msg::kDebug)
00598       << "New SimSnarlRecord " << simheader << endl;
00599     
00600     static VldTimeStamp rollbackTS;
00601     if ( fGetApplyRollback && mcGenTime < rollbackTS ) {
00602       rollbackTS = mcGenTime;
00603       // rollback the database for Ugli* and Bfld*
00604       DbiTableProxyRegistry& dbiCfg = DbiTableProxyRegistry::Instance();
00605       // something like: 
00606       //   dbiCfg.Set("Rollback:Ugli*         =  '2002-08-01'");
00607       dbiCfg.Set(Form("Rollback:UGLIDBI* = '%s'",rollbackTS.AsString("sql")));
00608       dbiCfg.Set(Form("Rollback:BFLDDBIPLANEMAP* = '%s'",rollbackTS.AsString("sql")));
00609       dbiCfg.Set("RollbackType:UGLIDBI* = 'CREATIONDATE'");
00610       dbiCfg.Set("RollbackType:BFLDDBIPLANEMAP* = 'CREATIONDATE'");
00611 
00612       MSG("Exodus",Msg::kInfo)
00613         << "RerootToTruthModule is setting dbiCfg to: " << endl;
00614       MsgService* msgserv = MsgService::Instance();
00615       if (msgserv->IsActive("Exodus",Msg::kInfo))
00616           dbiCfg.GetConfig().Print();
00617 
00618       dbiCfg.Update();
00619     }
00620 
00621     // give the RecSimSnarl to MOM to hold as a "fragment"
00622     record->GetTempTags().Set("stream","SimSnarl");  
00623     record->GetTempTags().Set("tree","SimSnarl");  
00624     record->GetTempTags().Set("file",RerootExodus::GetRerootFileName());
00625     record->GetTempTags().Set("index",RerootExodus::GetEventRecord());  
00626     mom->AdoptFragment(record);
00627   }
00628 
00629    // create/copy over other truth object
00630 
00631    // 2003-03-17 rwh: No longer copy REROOT_NeuVtx into record; doing
00632    // so is seriously flawed because there is a zoffset that depends
00633    // on the existence of gMINFast.  Thus one gets different values
00634    // for the REROOT_NeuVtx::Z() depending on whether one is running
00635    // from the original reroot file or a simsnarl'ed file.
00636 
00637    // Note: this next step is only legitimate if there is only ever 
00638    // *one* event in the snarl -- i.e. no overlay/pileup!
00639    REROOT_NeuKin *neukin_old = 
00640       (REROOT_NeuKin*)RerootExodus::GetNeuKinList()->At(0);
00641    if (neukin_old) record->AdoptComponent(new REROOT_NeuKin(*neukin_old));
00642 
00643    // Make a list of NeuKin objects
00644    // Keep the single object (done above) for backward compatibility
00645    record->AdoptComponent(make_neukin_list());
00646    record->AdoptComponent(make_fluxinfo_list());
00647    record->AdoptComponent(make_fluxwgt_list());
00648 
00649    // Add the StdHep list last for backward compatibility for those
00650    // that called SimSnarlRecord::FindComponent("TClonesArray") w/out
00651    // specifying the object name.  Set/Reset the PDG scheme used
00652    // for ions before/after conversion.
00653    UtilPDG::ionscheme_t old_ionscheme = UtilPDG::getDfltStdIonScheme();
00654    UtilPDG::setDfltStdIonScheme((UtilPDG::ionscheme_t)fIonScheme);
00655    // Make the list from what we have
00656    TClonesArray* stdhep_list = make_std_hep_list();
00657    ModifyStdHep(stdhep_list);
00658    if (fNewVtxFunction != "") RelocateStdHep(vldc,stdhep_list);
00659    record->AdoptComponent(stdhep_list);
00660    UtilPDG::setDfltStdIonScheme(old_ionscheme);
00661 
00662 #ifdef EVENT_KINEMATICS_PKGS
00663    // Add a list of NuEvtKin objects that link to StdHep
00664    record->AdoptComponent(make_NuEvtKin_list(stdhep_list));
00665 #endif
00666 
00667    // Add DigiScintHit's (if not disabled)
00668    if (fMakeDigiScintHitList)
00669      record->AdoptComponent(make_digi_scint_hit_list());
00670 
00671    // Create a DigiRerootInfo to hold miscellaneous information about
00672    // the event generation.
00673    DigiRerootInfo* drInfo = new DigiRerootInfo( RerootExodus::GetRunNo(),
00674                                                 RerootExodus::GetEventNo(),
00675                                                 RerootExodus::GetEventRecord(),
00676                                                 RerootExodus::BuildVldContext());
00677    record->AdoptComponent(drInfo);
00678 
00679    return JobCResult::kAOK; // All Ok
00680 }
00681 
00682 //......................................................................
00683 void RerootToTruthModule::ModifyStdHep(TClonesArray* stdhep_list)
00684 {
00685   // Fix up decay chain of StdHep list 
00686   // Trim it appropriately
00687 
00688   if (fFixCharmDKStatus) {
00689     // Fix cases where charm hadron generated particles (e.g. delta-rays)
00690     // before it decayed.  In this case the original was tagged ist=2,
00691     // the decay point was ist=1205 but the daughters had ist=205
00692     // rather than 1.  Be a bit careful in case we have charm->charm->stuff
00693     // (perhaps a D* -> D + gamma/pi0, but I haven't seen any of those
00694     // come out of NEUGEN in my study of a few files...)
00695 
00696     // decay products of charm that are ist=205 should be 1
00697     string charmDKDaughters1("fromcharm && ( ! ischarm ) && ist==205 ");
00698     UtilHepevt::modStatusStdHep(stdhep_list,charmDKDaughters1,1);
00699     // in case of charm(2)->charm(1205)->charm(205)->stuff, tag last charm 2
00700     string charmDKDaughters2("fromcharm && ischarm && ist==205 && da0 > -1  ");
00701     UtilHepevt::modStatusStdHep(stdhep_list,charmDKDaughters2,2);
00702   }
00703 
00704   if (fDropGeantEntries) {
00705     // drop entries added in GEANT propagation that aren't charm/tau related
00706     string dropGeant("(  (ist> 200 && ist< 231) || "
00707                      "   (ist> 300 && ist< 310) || "
00708                      "   (ist>1200 && ist<1231) || "
00709                      "   (ist>1300 && ist<1310) )  "
00710                      " && !ischarm && !istau ");
00711     UtilHepevt::modStatusStdHep(stdhep_list,dropGeant);
00712     UtilHepevt::squeezeStdHep(stdhep_list);
00713   }
00714   if (fDropTauDecayProducts) {
00715     // flag and remove tau decay point entry and decay daughters
00716     string dropTauDK("(istau && ist==1205) || (fromtau && ist==1)");
00717     UtilHepevt::modStatusStdHep(stdhep_list,dropTauDK);
00718     UtilHepevt::squeezeStdHep(stdhep_list);
00719     // reactivate tau as an active particle subject to decay IST=1
00720     UtilHepevt::modStatusStdHep(stdhep_list,"istau && ist==2 && da0==-1",1);
00721   }
00722   if (fDropCharmDecayProducts) {
00723      // flag and remove charm decay point entry and decay daughters
00724     string dropCharmDK("(ischarm && ist==1205) || (fromcharm && ist==1)");
00725     UtilHepevt::modStatusStdHep(stdhep_list,dropCharmDK);
00726     UtilHepevt::squeezeStdHep(stdhep_list);
00727     // reactivate tau as an active particle subject to decay IST=1
00728     UtilHepevt::modStatusStdHep(stdhep_list,"ischarm && ist==2 && da0==-1",1);
00729   }
00730   if (fDropStatus999) {
00731     UtilHepevt::modStatusStdHep(stdhep_list,"ist==999");
00732     UtilHepevt::squeezeStdHep(stdhep_list);
00733   }
00734   if (fStdHepEditString != "") {
00735     UtilHepevt::modStatusStdHep(stdhep_list,fStdHepEditString);
00736     UtilHepevt::squeezeStdHep(stdhep_list);
00737   }
00738 
00739 }
00740 
00741 //......................................................................
00742 void RerootToTruthModule::RelocateStdHep(const VldContext& vldc,
00743                                          TClonesArray* stdhep_list)
00744 {
00745   // Readjust the vertex, using fNewVtxFunction to do so
00746 
00747   MSG("Exodus",Msg::kDebug)
00748     << "RerootToTruthModule::RelocateStdHep() " << endl
00749     << "   vldc = " << vldc.AsString("c") << endl
00750     << "   fNewVtxFunction = \"" << fNewVtxFunction << "\"" << endl;
00751 
00752   int nhep = stdhep_list->GetLast()+1;
00753   if (nhep<1) return;  // nothing to do
00754 
00755   // Get the lead particle's vertex
00756   TParticle* leadpart = dynamic_cast<TParticle*>(stdhep_list->At(0));
00757   TLorentzVector vtx0;
00758   leadpart->ProductionVertex(vtx0);
00759 
00760   // Use user defined function to get new vertex
00761   // the signature needs to be of the form:
00762   //     void myfunction(VldContext*, TLorentzVector*)
00763   // where the function name is specified in fNewVtxFunction 
00764   // initialize the location of the new vertex with lead vtx
00765   // in case user has some use for that info ...
00766   TLorentzVector newvtx = vtx0;  
00767   TApplication* app = gROOT->GetApplication();
00768   std::string fmt = "";
00769   fmt += "VldContext* vldcp = (VldContext*)0x%lx; ";
00770   fmt += "TLorentzVector* newvtxp  = (TLorentzVector*)0x%lx; ";
00771   fmt += "%s(vldcp,newvtxp);";
00772   std::string cintline = 
00773     Form(fmt.c_str(),
00774          reinterpret_cast<unsigned long>(&vldc),
00775          reinterpret_cast<unsigned long>(&newvtx),
00776          fNewVtxFunction.c_str());
00777   MSG("Exodus",Msg::kDebug) << "ProcessLine \"" << cintline << "\"" << endl;
00778   app->ProcessLine(cintline.c_str());
00779   MSG("Exodus",Msg::kDebug) 
00780     << "got vertex " << endl
00781     << " vtx0   {" << vtx0[0] << "," << vtx0[1] 
00782     << ","<< vtx0[2] << "," << vtx0[3] << "}" << endl
00783     << " newvtx {" << newvtx[0] << "," << newvtx[1] 
00784     << ","<< newvtx[2] << "," << newvtx[3] << "}" << endl;
00785 
00786   // Apply new vertex to existing entries, by adding as a shift
00787   // to preserve relative offsets (e.g. charm/tau translation)
00788   // Avoid adding offsets to special cases such as non-particle lines, 
00789   // and positions relative to vtx for intranuke.
00790   TLorentzVector offset = newvtx - vtx0;
00791   for (int i = 0; i < nhep; ++i) {
00792     TParticle* apart = dynamic_cast<TParticle*>(stdhep_list->At(i));
00793     TLorentzVector partvtx;
00794     apart->ProductionVertex(partvtx);
00795     int istatus = apart->GetStatusCode();
00796 
00797     const Double_t r_mx = 0.001 * Munits::mm;
00798 
00799     bool status01 = istatus==0 || istatus==1 || istatus > 200;
00800     bool awayfrom0 = ( partvtx.Rho() > r_mx );
00801     if ( istatus!=999 && ( status01 || awayfrom0 ) ) {
00802       partvtx += offset;
00803       apart->SetProductionVertex(partvtx);
00804     } else {
00805       MSG("Exodus",Msg::kDebug) 
00806         << "no offset applied to ihep=" << i << endl
00807         << " status " << istatus 
00808         << " (test=" << (status01?"change":"x") << ") "
00809         << " nonzero vtx " << partvtx.Rho() << " vs. " << r_mx 
00810         << " (test=" << (awayfrom0?"change":"x") << ")" << endl;
00811     }
00812   }
00813 
00814 }
00815 
00816 //......................................................................
00817 JobCResult RerootToTruthModule::Ana(const MomNavigator *mom)
00818 {
00819    // Check out current RecSimSnarl Record from the current REROOT event
00820 
00821    const TObject* obj = 0;
00822 
00823    SimSnarlRecord* simrec = 
00824       dynamic_cast<SimSnarlRecord*>(mom->GetFragment("SimSnarlRecord"));
00825    if (!simrec) {
00826       MSG("Exodus", Msg::kWarning) << 
00827         "RerootToTruthModule::Ana no SimSnarlRecord" << endl;
00828       return JobCResult::kWarning;
00829    }
00830 
00831    // leading space
00832    cout << endl;
00833    // print the header
00834    simrec->GetSimSnarlHeader()->Print();
00835    cout << endl; // finish off header
00836 
00837    if (fAnaPrintStdHep) {
00838        obj = simrec->FindComponent("TClonesArray","StdHep");
00839        const TClonesArray* newstdhep = dynamic_cast<const TClonesArray*>(obj);
00840        if (!newstdhep) {
00841            MSG("Exodus", Msg::kWarning) << 
00842                "RerootToTruthModule::Ana no StdHep " << obj << endl;
00843            simrec->Print();
00844            return JobCResult::kWarning;
00845        }
00846        else print_tclones_stdhep(newstdhep);
00847    }
00848 
00849 
00850    // this only works out of the Reroot record (not the SimSnarlRecord)!
00851    if (fAnaPrintStdHepHead) print_tclones_stdhephead();
00852 
00853    if (fAnaPrintNeuKin) {
00854        obj = simrec->FindComponent("TClonesArray","NeuKinList");
00855        const TClonesArray* neukinlist = dynamic_cast<const TClonesArray*>(obj);
00856        if (neukinlist) {
00857            MSG("Exodus", Msg::kInfo) 
00858                << "RerootToTruthModule::Ana NeuKinList has "
00859                << neukinlist->GetEntries() << " entries " << endl;
00860            for (int ikin=0; ikin <= neukinlist->GetLast(); ikin++) {
00861                const REROOT_NeuKin* neukin = 
00862                    dynamic_cast<const REROOT_NeuKin*>(neukinlist->At(ikin));
00863                if ( neukin ) neukin->printEvent(cout);
00864                else cout << "NeuKinList entry " << ikin << " was empty!" << endl;
00865            }
00866        }
00867        else {
00868            MSG("Exodus", Msg::kWarning) << 
00869                "RerootToTruthModule::Ana no NeuKinList " << obj << endl;
00870            const REROOT_NeuKin* neukin = 
00871                dynamic_cast<const REROOT_NeuKin*>(simrec->FindComponent("REROOT_NeuKin"));
00872            if ( neukin ) neukin->printEvent(cout);
00873            else {
00874                simrec->Print();
00875                return JobCResult::kWarning;
00876            }
00877        }
00878    }
00879 
00880 #ifdef EVENT_KINEMATICS_PKGS
00881    if (fAnaPrintNuEvtKin) {
00882        obj = simrec->FindComponent("TClonesArray","NuEvtKinList");
00883        const TClonesArray* nuevtkinlist = dynamic_cast<const TClonesArray*>(obj);
00884        if (nuevtkinlist) {
00885            MSG("Exodus", Msg::kInfo) 
00886                << "RerootToTruthModule::Ana NuEvtKinList has "
00887                << nuevtkinlist->GetEntries() << " entries " << endl;
00888            for (int ik=0; ik <= nuevtkinlist->GetLast(); ik++) {
00889                const NuEvtKin* nuevtkin = 
00890                    dynamic_cast<const NuEvtKin*>(nuevtkinlist->At(ik));
00891                if ( nuevtkin ) nuevtkin->FormatToOStream(cout,"");  // "v" or ""
00892                else cout << "NuEvtKinList entry " << ik << " was empty!" << endl;
00893            }
00894        }
00895        else {
00896            MSG("Exodus", Msg::kWarning) << 
00897                "RerootToTruthModule::Ana no NuEvtKinList " << obj << endl;
00898            //simrec->Print();
00899            return JobCResult::kWarning;
00900        }
00901    }
00902 #endif
00903 
00904    if (fAnaPrintFluxInfo) {
00905        obj = simrec->FindComponent("TClonesArray","FluxInfoList");
00906        const TClonesArray* fluxinfolist = dynamic_cast<const TClonesArray*>(obj);
00907        if (fluxinfolist) {
00908            MSG("Exodus", Msg::kInfo) 
00909                << "RerootToTruthModule::Ana FluxInfoList has "
00910                << fluxinfolist->GetEntries() << " entries " << endl;
00911            for (int ikin=0; ikin <= fluxinfolist->GetLast(); ikin++) {
00912                const REROOT_FluxInfo* fluxinfo = 
00913                    dynamic_cast<const REROOT_FluxInfo*>(fluxinfolist->At(ikin));
00914                if ( fluxinfo ) fluxinfo->printEvent(cout);
00915                else cout << "FluxInfoList entry " << ikin << " was empty!" << endl;
00916            }
00917        }
00918        else {
00919            MSG("Exodus", Msg::kWarning) << 
00920                "RerootToTruthModule::Ana no FluxInfoList " << obj << endl;
00921        }
00922    }
00923 
00924    if (fAnaPrintFluxWgt) {
00925        obj = simrec->FindComponent("TClonesArray","FluxWgtList");
00926        const TClonesArray* fluxwgtlist = dynamic_cast<const TClonesArray*>(obj);
00927        if (fluxwgtlist) {
00928            MSG("Exodus", Msg::kInfo) 
00929                << "RerootToTruthModule::Ana FluxWgtList has "
00930                << fluxwgtlist->GetEntries() << " entries " << endl;
00931            for (int ikin=0; ikin <= fluxwgtlist->GetLast(); ikin++) {
00932                const REROOT_FluxWgt* fluxwgt = 
00933                    dynamic_cast<const REROOT_FluxWgt*>(fluxwgtlist->At(ikin));
00934                if ( fluxwgt ) fluxwgt->printEvent(cout);
00935                else cout << "FluxWgtList entry " << ikin << " was empty!" << endl;
00936            }
00937        }
00938        else {
00939            MSG("Exodus", Msg::kWarning) << 
00940                "RerootToTruthModule::Ana no FluxWgtList " << obj << endl;
00941        }
00942    }
00943 
00944    if (fAnaPrintScintHit) {
00945        obj = simrec->FindComponent("TClonesArray","DigiScintHits");
00946        const TClonesArray* scinthitlist = dynamic_cast<const TClonesArray*>(obj);
00947        if (scinthitlist) {
00948            MSG("Exodus", Msg::kInfo) 
00949                << "RerootToTruthModule::Ana ScintHitList has "
00950                << scinthitlist->GetEntries() << " entries " << endl;
00951            for (int ihit=0; ihit <= scinthitlist->GetLast(); ihit++) {
00952                const DigiScintHit* scinthit = 
00953                    dynamic_cast<const DigiScintHit*>(scinthitlist->At(ihit));
00954                if ( scinthit ) scinthit->Print();
00955                else cout << "ScintHitList entry " << ihit << " was empty!" << endl;
00956            }
00957        }
00958        else {
00959            MSG("Exodus", Msg::kWarning) << 
00960                "RerootToTruthModule::Ana no ScintHitList " << obj << endl;
00961        }
00962    }
00963 
00964 
00965    return JobCResult::kPassed; // All Ok
00966 }
00967 //......................................................................
00968 
00969 const Registry& RerootToTruthModule::DefaultConfig() const
00970 {
00971 //======================================================================
00972 // Create a registry which holds the default configuration and return it
00973 //======================================================================
00974   static Registry r;
00975  
00976   // Set name of config
00977   std::string name = this->GetName();
00978   name += ".config.default";
00979   r.SetName(name.c_str());
00980 
00981   // Set values of config
00982   r.UnLockValues();
00983   r.Set("GetApplyRollback",    1);
00984   r.Set("AnaPrintStdHep",      1);
00985   r.Set("AnaPrintStdHepHead",  0);
00986   r.Set("AnaPrintNeuKin",      1);
00987   r.Set("AnaPrintNuEvtKin",    1);
00988   r.Set("AnaPrintFluxInfo",    1);
00989   r.Set("AnaPrintFluxWgt",     1);
00990   r.Set("AnaPrintScintHit",    0);
00991 
00992   r.Set("VldSimFlag","kUnknown");
00993 
00994   r.Set("ConvertIonPDG",(int)UtilPDG::kIonUnchanged);
00995   r.Set("FixCharmDKStatus",       1);  // do *this* by default
00996   r.Set("DropGeantEntries",       0);
00997   r.Set("DropTauDecayProducts",   0);
00998   r.Set("DropCharmDecayProducts", 0);
00999   r.Set("DropStatus999",          0);
01000   r.Set("StdHepEditString",      "");
01001 
01002   r.Set("MakeDigiScintHitList",   fgGlobalMakeDigiScintList);
01003 
01004   r.LockValues();
01005 
01006   return r;
01007 }
01008 
01009 //......................................................................
01010 
01011 void RerootToTruthModule::Config(const Registry& r)
01012 {
01013 //======================================================================
01014 // Configure the module given the registry r
01015 //======================================================================
01016 //  char   tmpb;
01017 //  double tmpd;
01018   int    tmpi;
01019 
01020 
01021   if (r.Get("GetApplyRollback",    tmpi)) { fGetApplyRollback   = tmpi; }
01022   if (r.Get("AnaPrintStdHep",      tmpi)) { fAnaPrintStdHep     = tmpi; }
01023   if (r.Get("AnaPrintStdHepHead",  tmpi)) { fAnaPrintStdHepHead = tmpi; }
01024   if (r.Get("AnaPrintNeuKin",      tmpi)) { fAnaPrintNeuKin     = tmpi; }
01025   if (r.Get("AnaPrintNuEvtKin",    tmpi)) { fAnaPrintNuEvtKin   = tmpi; }
01026   if (r.Get("AnaPrintFluxInfo",    tmpi)) { fAnaPrintFluxInfo   = tmpi; }
01027   if (r.Get("AnaPrintFluxWgt",     tmpi)) { fAnaPrintFluxWgt    = tmpi; }
01028   if (r.Get("AnaPrintScintHit",    tmpi)) { fAnaPrintScintHit   = tmpi; }
01029 
01030   if (r.Get("ConvertIonPDG",tmpi)) {
01031     if (fIonScheme != tmpi) {
01032       MSG("Exodus",Msg::kInfo) 
01033         << GetName() << " will convert ions using the "
01034         << UtilPDG::ionSchemeName((UtilPDG::ionscheme_t)tmpi)
01035         << " scheme." << endl;
01036     }
01037     fIonScheme = tmpi;
01038   }
01039 
01040   if (r.Get("FixCharmDKStatus",      tmpi)) { fFixCharmDKStatus       = tmpi; }
01041   if (r.Get("DropGeantEntries",      tmpi)) { fDropGeantEntries       = tmpi; }
01042   if (r.Get("DropTauDecayProducts",  tmpi)) { fDropTauDecayProducts   = tmpi; }
01043   if (r.Get("DropCharmDecayProducts",tmpi)) { fDropCharmDecayProducts = tmpi; }
01044   if (r.Get("DropStatus999",         tmpi)) { fDropStatus999          = tmpi; }
01045   fStdHepEditString = r.GetCharString("StdHepEditString");
01046 
01047   if (r.Get("MakeDigiScintHitList",  tmpi)) { fMakeDigiScintHitList   = tmpi; }
01048 
01049   MSG("Exodus",Msg::kInfo)
01050     << GetName() << " StdHep list editing,"
01051     << " FixCharmDKStatus " << (fFixCharmDKStatus?"yes":"no")
01052     <<", Drop " << endl
01053     << " GeantEntries="       << (fDropGeantEntries?"yes":"no")
01054     << " TauDecayProducts="   << (fDropTauDecayProducts?"yes":"no")
01055     << " CharmDecayProducts=" << (fDropCharmDecayProducts?"yes":"no")
01056     << " Status999="          << (fDropStatus999?"yes":"no")
01057     << endl;
01058   if (fStdHepEditString != "") {
01059     MSG("Exodus",Msg::kInfo)
01060       << " StdHepEditString: \"" << fStdHepEditString << "\""
01061       << endl;
01062   }
01063   MSG("Exodus",Msg::kInfo)
01064     << GetName() << " MakeDigiScintHitList = " 
01065     << (fMakeDigiScintHitList?"yes":"no") << endl;
01066 
01067   const char* simflg_string = r.GetCharString("VldSimFlag");
01068   SimFlag::SimFlag_t simflg = SimFlag::StringToEnum(simflg_string);
01069   if ( simflg == SimFlag::kUnknown ) simflg = RerootExodus::GetVldSimFlag();
01070   if ( simflg != RerootExodus::GetVldSimFlag() ) {
01071     MSG("Exodus",Msg::kInfo)
01072       << GetName() << " Changing SimFlag from "
01073       << SimFlag::AsString(RerootExodus::GetVldSimFlag()) << " to "
01074       << SimFlag::AsString(simflg) << " in SimSnarlHeader creation."
01075       << endl;
01076     RerootExodus::SetVldSimFlag(simflg);
01077   } else {
01078     MSG("Exodus",Msg::kInfo)
01079       << GetName() << " Using SimFlag "
01080       << SimFlag::AsString(simflg) << " in SimSnarlHeader creation."
01081       << endl;
01082   }
01083 
01084 }
01085 
01086 //......................................................................