RerootToTruthModule.cxx File Reference

#include <cassert>
#include <RerootExodus/RerootToTruthModule.h>
#include <RerootExodus/RerootExodus.h>
#include <MinosObjectMap/MomNavigator.h>
#include <MessageService/MsgService.h>
#include <JobControl/JobCModuleRegistry.h>
#include <Conventions/Mphysical.h>
#include <REROOT_Classes/REROOT_StdHep.h>
#include <REROOT_Classes/REROOT_StdHepHead.h>
#include <REROOT_Classes/REROOT_NeuKin.h>
#include <REROOT_Classes/REROOT_FluxInfo.h>
#include <REROOT_Classes/REROOT_FluxWgt.h>
#include <Record/SimSnarlRecord.h>
#include <Record/SimSnarlHeader.h>
#include <UgliGeometry/UgliGeomHandle.h>
#include <UgliGeometry/UgliStripHandle.h>
#include <Digitization/DigiScintHit.h>
#include <Digitization/DigiRerootInfo.h>
#include <DatabaseInterface/DbiTableProxyRegistry.h>
#include <Util/LoadMinosPDG.h>
#include <Util/UtilPDG.h>
#include <Util/UtilHepevt.h>
#include <TMath.h>
#include <TString.h>
#include <TSystem.h>
#include <TDatabasePDG.h>
#include <TParticle.h>
#include <TLorentzVector.h>
#include <TROOT.h>
#include <TApplication.h>

Go to the source code of this file.

Functions

 CVSID ("$Id: RerootToTruthModule.cxx,v 1.51 2008/05/16 15:50:29 rhatcher Exp $")
 JOBMODULE (RerootToTruthModule,"RerootToTruthModule","Builds RawRecord with RawDigitDataBlock")
static TClonesArray * make_std_hep_list ()
static TClonesArray * make_neukin_list ()
static TClonesArray * make_fluxinfo_list ()
static TClonesArray * make_fluxwgt_list ()
static double flshit_path_length (REROOT_FLSHit *h)
static double flshit_track_energy (REROOT_FLSHit *h)
static void flshit_massage_local (PlexStripEndId seid, float x, float y, float z, float &xout, float &yout, float &zout)
static TClonesArray * make_digi_scint_hit_list ()
static void print_tclones_stdhep (const TClonesArray *newstdhep)
static void print_tclones_stdhephead ()

Variables

static std::string fgGlobalNewVtxFunction = ""
static int fgGlobalMakeDigiScintList = 1

Function Documentation

CVSID ( "$Id: RerootToTruthModule.  cxx,
v 1.51 2008/05/16 15:50:29 rhatcher Exp $"   
)
static void flshit_massage_local ( PlexStripEndId  seid,
float  x,
float  y,
float  z,
float &  xout,
float &  yout,
float &  zout 
) [static]

Definition at line 335 of file RerootToTruthModule.cxx.

References RerootExodus::BuildVldContext(), Munits::cm, PlexPlaneId::GetNumStrips(), PlexStripEndId::GetStrip(), UgliGeomHandle::GetStripHandle(), UgliStripHandle::GlobalToLocal(), UgliGeomHandle::IsValid(), Msg::kWarning, MSG, and RerootExodus::SEIdLocalToGlobal().

Referenced by make_digi_scint_hit_list().

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 }

static double flshit_path_length ( REROOT_FLSHit h  )  [static]

Definition at line 298 of file RerootToTruthModule.cxx.

References Munits::cm, REROOT_FLSHit::XBegin(), REROOT_FLSHit::XEnd(), REROOT_FLSHit::YBegin(), REROOT_FLSHit::YEnd(), REROOT_FLSHit::ZBegin(), and REROOT_FLSHit::ZEnd().

Referenced by make_digi_scint_hit_list().

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 }

static double flshit_track_energy ( REROOT_FLSHit h  )  [static]

Definition at line 305 of file RerootToTruthModule.cxx.

References Munits::GeV, REROOT_FLSHit::IPDG(), Msg::kWarning, Munits::m, MSG, and REROOT_FLSHit::Ptot().

Referenced by make_digi_scint_hit_list().

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 }

JOBMODULE ( RerootToTruthModule  ,
"RerootToTruthModule"  ,
"Builds RawRecord with RawDigitDataBlock  
)
static TClonesArray* make_digi_scint_hit_list (  )  [static]

Definition at line 378 of file RerootToTruthModule.cxx.

References Munits::c_light, Munits::cm, MuELoss::e, REROOT_FLSHit::ELoss(), flshit_massage_local(), flshit_path_length(), flshit_track_energy(), RerootExodus::GetFLSHitList(), REROOT_FLSHit::ICell(), REROOT_FLSHit::IExtr(), REROOT_FLSHit::IPDG(), REROOT_FLSHit::IPln(), REROOT_FLSHit::ITrack(), Msg::kWarning, StripEnd::kWest, MSG, RerootExodus::PECAB2SEId(), RerootExodus::PlaneName(), PlexStripEndId::SetEnd(), size, REROOT_FLSHit::TOFG(), REROOT_FLSHit::XBegin(), REROOT_FLSHit::XEnd(), REROOT_FLSHit::YBegin(), REROOT_FLSHit::YEnd(), REROOT_FLSHit::ZBegin(), and REROOT_FLSHit::ZEnd().

Referenced by RerootToTruthModule::Get().

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 }

static TClonesArray* make_fluxinfo_list (  )  [static]

Definition at line 179 of file RerootToTruthModule.cxx.

References RerootExodus::GetFluxInfoList().

Referenced by RerootToTruthModule::Get().

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 }

static TClonesArray* make_fluxwgt_list (  )  [static]

Definition at line 199 of file RerootToTruthModule.cxx.

References RerootExodus::GetFluxWgtList().

Referenced by RerootToTruthModule::Get().

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 }

static TClonesArray* make_neukin_list (  )  [static]

Definition at line 159 of file RerootToTruthModule.cxx.

References RerootExodus::GetNeuKinList().

Referenced by RerootToTruthModule::Get().

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 }

static TClonesArray* make_std_hep_list (  )  [static]

Definition at line 117 of file RerootToTruthModule.cxx.

References Mphysical::c_light, REROOT_StdHep::E(), REROOT_StdHep::FirstChild(), REROOT_StdHep::FirstParent(), RerootExodus::GetStdHepList(), Munits::GeV, REROOT_StdHep::IdHEP(), REROOT_StdHep::IstHEP(), REROOT_StdHep::LastChild(), REROOT_StdHep::LastParent(), Munits::mm, REROOT_StdHep::Px(), REROOT_StdHep::Py(), REROOT_StdHep::Pz(), UtilPDG::stdIonNumber(), REROOT_StdHep::Tprod(), REROOT_StdHep::Xmm(), REROOT_StdHep::Ymm(), and REROOT_StdHep::Zmm().

Referenced by RerootToTruthModule::Get().

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 }

static void print_tclones_stdhep ( const TClonesArray *  newstdhep  )  [static]

Definition at line 451 of file RerootToTruthModule.cxx.

References Msg::kInfo, and MSG.

Referenced by RerootToTruthModule::Ana().

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 }

static void print_tclones_stdhephead (  )  [static]

Definition at line 525 of file RerootToTruthModule.cxx.

References RerootExodus::GetStdHepHeadList(), REROOT_StdHepHead::ID(), REROOT_StdHepHead::NevHEP(), and REROOT_StdHepHead::NHEP().

Referenced by RerootToTruthModule::Ana().

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 }


Variable Documentation

int fgGlobalMakeDigiScintList = 1 [static]
std::string fgGlobalNewVtxFunction = "" [static]

Generated on 15 Nov 2018 for loon by  doxygen 1.6.1