loon: RerootToTruthModule.cxx Source File

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