MCMonitorCosmicHistograms.cxx

Go to the documentation of this file.

//
// Histograms used by MCMonitor (includes useful ancillary functions)
//
// Kregg Arms (arms@physics.umn.edu)
#include <iostream>
#include <sstream>
#include "TClonesArray.h"
#include "TDatabasePDG.h"
#include "TParticlePDG.h"

#include "MCMonitorCosmicHistograms.h"
#include "MessageService/MsgService.h"
#include "Conventions/Munits.h"
#include "HistMan/HistMan.h"
#include "MCNtuple/NtpMCRecord.h"
#include "MCNtuple/NtpMCDigiScintHit.h"
#include "MCNtuple/NtpMCStdHep.h"
#include "MCNtuple/NtpMCTruth.h"
#include "CandNtupleSR/NtpSRRecord.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "TruthHelperNtuple/NtpTHRecord.h"
#include "TruthHelperNtuple/NtpTHTrack.h"
#include "StandardNtuple/NtpStRecord.h"
#include "Plex/PlexPlaneId.h"

CVSID("$Id: MCMonitorCosmicHistograms.cxx,v 1.1 2007/10/10 03:11:22 schubert Exp $");

//_______________________________________________________________________
MCMonitorCosmicHistograms::MCMonitorCosmicHistograms(std::string name)
: fName(name),fDetector(Detector::kUnknown),fShieldFlag(kAll),fParticleId(0) {
  // Default constructor

}

//_______________________________________________________________________
MCMonitorCosmicHistograms::~MCMonitorCosmicHistograms() {
  // Destructor

}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillHistograms(const NtpStRecord* strec,
                                               const NtpMCRecord* mcrec,
                                               const NtpSRRecord* srrec,
                                               const NtpTHRecord* threc) {
  // Public method to fill all histograms. Histograms created on first call.

  if (!strec && !mcrec && !srrec && !threc) return;
  
  // Set detector type
  const VldContext* vldc 
           = (strec)? strec->GetVldContext() : mcrec->GetVldContext();
  fDetector = vldc->GetDetector();
  
  std::string basepath = fName;
  FillHistsTrue(basepath,strec,mcrec);
  FillHistsReco(basepath,strec,srrec);
  FillHistsReTr(basepath,strec,mcrec,srrec,threc);

}
  
//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillHistsTrue(std::string basepath,
                                              const NtpStRecord* strec,
                                              const NtpMCRecord* mcrec) {
  // Private method to fill all ../True/.. histograms

  if (!strec && !mcrec) return;
  
  std::string truepath = basepath + "/True";
  
  const NtpMCSummary& mchdr = (strec)? strec->mchdr : mcrec->mchdr;
  
  // mchdr
  FillTrueNumEvt(truepath,mchdr.nmc);
  FillTrueNumStdHep(truepath,mchdr.nstdhep);
  FillTrueNumDigiHit(truepath,mchdr.ndigihit);

  const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
  const int npid = 5;
  int pidarray[npid] = {13,-13,11,-11,2212}; // mu-,mu+,e-,e+,p
  
  // digihit
  const TClonesArray* digihitarray 
                       = (strec) ? strec->digihit : mcrec->digihit;
  if ( digihitarray ) {
    // Fill first without shield hits
    fShieldFlag = kNoShield;
    fParticleId = 0;
    FillHistsTrueDigiHit(truepath+"/NoShield",digihitarray,
                        " (VS Excl)");
    // Now fill for each particle type of interest
    for ( int ipid = 0; ipid < npid; ipid++ ) {
      fParticleId = pidarray[ipid]; // mu+
      std::string partname = "???";
      if ( dbpdg.GetParticle(fParticleId) ) partname
                              = dbpdg.GetParticle(fParticleId)->GetName();
      std::string partpath = truepath+"/NoShield/"+partname;
      std::string clause = " (VS Excl) (" + partname + ")";
      FillHistsTrueDigiHit(partpath,digihitarray,clause);
    }
    
    // Fill next with shield hits only
    fParticleId = 0;
    fShieldFlag = kShieldOnly;
    FillHistsTrueDigiHit(truepath+"/ShieldOnly",digihitarray,
                                               " (VS Only)");
    fShieldFlag = kAll;
  }
  
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillHistsReco(std::string basepath,
                                              const NtpStRecord* strec,
                                              const NtpSRRecord* srrec) {
  // Private method to fill all ../Reco/.. histograms

  if (!strec && !srrec) return;
  
  std::string recopath = basepath + "/Reco";
  
  // strip
  const TClonesArray* stparray = (strec) ? strec->stp: srrec->stp;
  FillHistsRecoStp(recopath,stparray);
  
  // shower
  const TClonesArray* shwarray = (strec) ? strec->shw : srrec->shw;
  FillHistsRecoShw(recopath,shwarray);

  // track
  const TClonesArray* trkarray = (strec) ? strec->trk : srrec->trk;
  FillHistsRecoTrk(recopath,trkarray);

  // event
  //const TClonesArray* evtarray = (strec) ? strec->evt : srrec->evt;
  //FillHistsRecoEvt(recopath,evtarray);
  
  
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillHistsReTr(std::string basepath,
                                              const NtpStRecord* strec,
                                              const NtpMCRecord* mcrec,
                                              const NtpSRRecord* srrec,
                                              const NtpTHRecord* threc) {
  // Private method to fill all ../ReTr/.. histograms

  if (!strec && (!mcrec || !srrec || !threc) ) return;
  
  std::string retrpath = basepath + "/ReTr";

  // tracks
  const TClonesArray* stdheparray = (strec) ? strec->stdhep : mcrec->stdhep;
  const TClonesArray* trkarray = (strec) ? strec->trk : srrec->trk;
  const TClonesArray* thtrkarray = (strec) ? strec->thtrk : threc->thtrk;
  FillHistsReTrTrk(retrpath,stdheparray,trkarray,thtrkarray);
  
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillHistsReTrTrk(std::string retrpath,
                                              const TClonesArray* stdheparray,
                                              const TClonesArray* trkarray,
                                              const TClonesArray* thtrkarray) {
  // Private method to fill all ../ReTr/(PartType)/Trk histograms
  
  if ( !stdheparray || !trkarray || !thtrkarray ) return;
  
  for (int itrk = 0; itrk < thtrkarray->GetEntriesFast(); itrk++) {
    NtpTHTrack* thtrk = dynamic_cast<NtpTHTrack*>(thtrkarray->At(itrk));
    Int_t istdhep = thtrk->trkstdhep;
    if ( istdhep < 0 ) continue; // truth helper failed to find match if true
    NtpMCStdHep* stdhep=dynamic_cast<NtpMCStdHep*>(stdheparray->At(istdhep));
    std::string partname = stdhep->ParticleName();
    std::string partpath = retrpath + "/" + partname;
    NtpSRTrack*  srtrk = dynamic_cast<NtpSRTrack*>(trkarray->At(itrk));
    std::string clause = " (" + partname + ")";
    FillReTrDelQPCurve(partpath,*stdhep,*srtrk,clause);
  }

}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillReTrDelQPCurve(std::string partpath,
                            const NtpMCStdHep& stdhep,const NtpSRTrack& trk,
                            std::string clause) {
  // Private method to fill ../DelQPCurve histogram

  HistMan hm(partpath.c_str());
  TH1F* hist1 = hm.Get<TH1F>("DelQPCurve");

  if ( !hist1 ) {
    hist1 = new TH1F("DelQPCurve","[p/q - p/q(curvature)]",40,-100.,100.);
    std::string xtitle = "[p/q(true) - p/q(curvature)]"+clause;
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/5. GeV/c");
    hm.Adopt("",hist1);
  }
  
  double trkqp = (TMath::Abs(trk.momentum.qp) > 0.) ? (1/trk.momentum.qp) : 0;
  double stdhepqp = (stdhep.Charge() != 0 ) ? stdhep.P()/stdhep.Charge() : 0;
  hist1 -> Fill(stdhepqp - trkqp);
                                              
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillHistsRecoStp(std::string recopath,
                                         const TClonesArray* stparray) {
  // Private method to fill ../Stp* histograms
  
  if ( !stparray ) return;

  HistMan hm(recopath.c_str());

  // For performance reasons, pre-read stp hists
  TH1F* eastphpe = 0;
  TH1F* westphpe = 0;
  FillRecoStpPerSnarl(recopath,stparray->GetEntriesFast());
  
  for ( int istp = 0; istp < stparray->GetEntriesFast(); istp++ ) {
    const NtpSRStrip& strip 
             = *(dynamic_cast<NtpSRStrip*>(stparray->At(istp)));

    if ( !eastphpe ) eastphpe = hm.Get<TH1F>("StpEastPHPE");
    if ( !westphpe ) westphpe = hm.Get<TH1F>("StpWestPHPE");
    
    // Determine if this is shield associated strip or not
    PlexPlaneId plexId(fDetector,strip.plane);
    bool isShield = plexId.IsVetoShield();
    if ( isShield ) continue;

    FillRecoStpEastPHPE(recopath,strip,eastphpe," (VS Excl) ");
    FillRecoStpWestPHPE(recopath,strip,westphpe," (VS Excl) ");
 
  }

  return;
  
}


//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillHistsRecoTrk(std::string recopath,
                                         const TClonesArray* trkarray) {
  // Private method to fill all ../Reco/Trk* histograms
  
  if ( !trkarray ) return;

  HistMan hm(recopath.c_str());
  FillRecoTrkPerSnarl(recopath,trkarray->GetEntriesFast());
  
  for ( int itrk = 0; itrk < trkarray->GetEntriesFast(); itrk++ ) {
    const NtpSRTrack& track 
             = *(dynamic_cast<NtpSRTrack*>(trkarray->At(itrk)));

    FillRecoTrkVtxDCosX(recopath,track);
    FillRecoTrkVtxDCosY(recopath,track);
    FillRecoTrkVtxDCosZ(recopath,track);
    FillRecoTrkRange(recopath,track);
    
  }
  
  return;
  
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillRecoStpEastPHPE(std::string recopath,
                                                    const NtpSRStrip& stp,
                                                    TH1F* hist1,
                                                    std::string clause) {
  // Private method to fill ../StpEastPHPE histogram

  if ( !hist1 ) {
    HistMan hm(recopath.c_str());
    hist1 = new TH1F("StpEastPHPE","Strip East PH",100,0.,2000.);
    std::string xtitle = "Strip East PH(PE)" + clause;
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/20 PE");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(stp.ph0.sigcor);
                                              
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillRecoStpWestPHPE(std::string recopath,
                                                    const NtpSRStrip& stp,
                                                    TH1F* hist1,
                                                    std::string clause) {
  // Private method to fill ../StpWestPHPE histogram

  if ( !hist1 ) {
    HistMan hm(recopath.c_str());
    hist1 = new TH1F("StpWestPHPE","Strip West PH",100,0.,2000.);
    std::string xtitle = "Strip West PH(PE)" + clause;
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/20 PE");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(stp.ph1.sigcor);
                                              
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillRecoTrkVtxDCosX(std::string recopath,
                                                    const NtpSRTrack& trk,
                                                    std::string particlename) {
  // Private method to fill ../TrkVtxDCosX histogram

  HistMan hm(recopath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("TrkVtxDCosX");
  if ( !hist1 ) {
    hist1 = new TH1F("TrkVtxDCosX","reco trk vtx dcos(x)",40,-1.,1.);
    std::string xtitle = particlename + " Reco Trk Vtx DCos(X)";
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/0.05");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(trk.vtx.dcosx);
                                              
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillRecoTrkVtxDCosY(std::string recopath,
                                                    const NtpSRTrack& trk,
                                                    std::string particlename) {
  // Private method to fill all ../TrkVtxDCosY histogram

  HistMan hm(recopath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("TrkVtxDCosY");
  if ( !hist1 ) {
    hist1 = new TH1F("TrkVtxDCosY","reco trk vtx dcos(y)",40,-1.,1.);
    std::string xtitle = particlename + " Reco Trk Vtx DCos(Y)";
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/0.05");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(trk.vtx.dcosy);
                                              
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillRecoTrkVtxDCosZ(std::string recopath,
                                                    const NtpSRTrack& trk,
                                                    std::string particlename) {
  // Private method to fill all ../TrkVtxDCosZ histograms

  HistMan hm(recopath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("TrkVtxDCosZ");
  if ( !hist1 ) {
    hist1 = new TH1F("TrkVtxDCosZ","reco trk vtx dcos(z)",40,-1.,1.);
    std::string xtitle = particlename + " Reco Trk Vtx DCos(Z)";
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/0.05");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(trk.vtx.dcosz);
                                              
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillRecoTrkRange(std::string recopath,
                                                 const NtpSRTrack& trk,
                                                 std::string particlename) {
  // Private method to fill ../TrkRange histogram

  HistMan hm(recopath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("TrkRange");
  if ( !hist1 ) {
    hist1 = new TH1F("TrkRange","reco trk range",64,0.,32.);
    std::string xtitle = particlename + " Reco Trk Range (m)";
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/0.5 m");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(trk.ds);
                                              
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillRecoTrkPerSnarl(std::string recopath,
                                                    unsigned int ntrack) {
  // Private method to fill ../TrkPerSnarl histogram

  HistMan hm(recopath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("TrkPerSnarl");
  if ( !hist1 ) {
    hist1 = new TH1F("TrkPerSnarl","track per snarl",21,-0.5,20.5);
    hist1->SetXTitle("tracks/snarl");
    hist1->SetYTitle("Entries/1");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(ntrack);
                                              
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillRecoShwPerSnarl(std::string recopath,
                                                    unsigned int nshower) {
  // Private method to fill ../ShwPerSnarl histogram

  HistMan hm(recopath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("ShwPerSnarl");
  if ( !hist1 ) {
    hist1 = new TH1F("ShwPerSnarl","shower per snarl",21,-0.5,20.5);
    hist1->SetXTitle("showers/snarl");
    hist1->SetYTitle("Entries/1");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(nshower);
                                              
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillRecoStpPerSnarl(std::string recopath,
                                                    unsigned int nstrip) {
  // Private method to fill ../StpPerSnarl histogram

  HistMan hm(recopath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("StpPerSnarl");
  if ( !hist1 ) {
    hist1 = new TH1F("StpPerSnarl","strip per snarl",100,-0.5,999.5);
    hist1->SetXTitle("strips/snarl");
    hist1->SetYTitle("Entries/1");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(nstrip);
                                              
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillHistsRecoShw(std::string recopath,
                                       const TClonesArray* shwarray) {
  // Private method to fill all ../Reco/Shw* histograms
  
  if ( !shwarray ) return;

  HistMan hm(recopath.c_str());

  FillRecoShwPerSnarl(recopath,shwarray->GetEntriesFast());
  
  for ( int ishw = 0; ishw < shwarray->GetEntriesFast(); ishw++ ) {
    const NtpSRShower& shower 
             = *(dynamic_cast<NtpSRShower*>(shwarray->At(ishw)));

    FillRecoShwNumPlanes(recopath,shower);
    FillRecoShwEMGeV(recopath,shower);
    
  }
  
  return;
  
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillRecoShwNumPlanes(std::string recopath,
                                              const NtpSRShower& shw) {
  // Private method to fill all ../ShwNumPlanes histograms

  HistMan hm(recopath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("ShwNumPlanes");
  if ( !hist1 ) {
    hist1 = new TH1F("ShwNumPlanes","Shw Num Planes",25,0.,50.);
    hist1->SetXTitle("Shwr Num Planes");
    hist1->SetYTitle("Entries/2");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(shw.plane.n);
  
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillRecoShwEMGeV(std::string recopath,
                                                 const NtpSRShower& shw) {

  // Private method to fill all ../ShwEMGeV histograms

  HistMan hm(recopath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("ShwEMGeV");
  if ( !hist1 ) {
    hist1 = new TH1F("ShwEMGeV","Shw EM",40,0.,20.);
    hist1->SetXTitle("Shwr EM (GeV)");
    hist1->SetYTitle("Entries/0.5 GeV");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(shw.shwph.EMgev);
  
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueNumEvt(std::string truepath,
                                               int nmc) {
  // Private method to fill ../True/NumEvtPerSnarl 
  
  HistMan hm(truepath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("NumEvtPerSnarl");
  if ( !hist1 ) {
    hist1 = new TH1F("NumEvtPerSnarl","True Num Evt/Snarl",60,0.,60.);
    hist1->SetXTitle("True Num Evt/Snarl");
    hist1->SetYTitle("Entries/1");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(nmc);
  
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueNumStdHep(std::string truepath,
                                                  int nstdhep) {
  // Private method to fill ../True/NumStdHepPerSnarl 
  
  HistMan hm(truepath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("NumStdHepPerSnarl");
  if ( !hist1 ) {
    hist1 = new TH1F("NumStdHepPerSnarl","Num StdHep/Snarl",50,0.,100.);
    hist1->SetXTitle("Num StdHep/Snarl");
    hist1->SetYTitle("Entries/2");
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(nstdhep);
  
}
  
//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueNumDigiHit(std::string truepath,
                                      int ndigihit, std::string clause) {
  // Private method to fill ../True/NumDigiHitPerSnarl 
  
  HistMan hm(truepath.c_str());

  TH1F* hist1 = hm.Get<TH1F>("NumDigiHitPerSnarl");
  if ( !hist1 ) {
    Float_t maxrange = 10000;
    Int_t nbin = 50;
    if ( fParticleId == -13 || fParticleId == 13 || fParticleId == -11 ) 
                                             { maxrange=250; nbin=50; }
    //if ( fParticleId == 2212 ) { maxrange = 10; nbin = 40; }
    
    hist1 = new TH1F("NumDigiHitPerSnarl","Num DigiHit/Snarl",nbin,
                                                      0.,maxrange);
    std::string xtitle = "Num DigiHits/Snarl" + clause;
    hist1->SetXTitle(xtitle.c_str());
    Float_t binwidth = maxrange/nbin;
    std::ostringstream oss;
    oss << std::dec << binwidth;
    std::string ytitle = "Entries/"+oss.str();
    hist1->SetYTitle(ytitle.c_str());
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(ndigihit);
  
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillHistsTrueDigiHit(std::string truepath,
                                       const TClonesArray* digihitarray,
                                       std::string clause) {
  // Private method to fill all ../True/DigiHit/.. histograms
  
  if ( !digihitarray ) return;
  std::string digihitpath = truepath + "/DigiHit";

  HistMan hm(digihitpath.c_str());
  // For performance reasons, pre-read digihit hists
  TH1F* dshist1 = 0;
  TH1F* dsbyephist1 = 0;
  TProfile* dsbyepvsdsprof = 0;
  TH1F* dehist1 = 0;
  TH1F* dedshist1 = 0;
  TH1F* dedsbyephist1 = 0;
  TProfile* dedsvspeprof = 0;
  TProfile* dedsbyepvspeprof = 0;
  
  int ndigihit = 0;
  for ( int ihit = 0; ihit < digihitarray->GetEntriesFast(); ihit++ ) {
    const NtpMCDigiScintHit& digihit 
             = *(dynamic_cast<NtpMCDigiScintHit*>(digihitarray->At(ihit)));

    // Determine if this is shield associated hit or not
    PlexPlaneId plexId(fDetector,digihit.plane);
    bool isShield = plexId.IsVetoShield();
    if ( isShield && fShieldFlag == kNoShield ) continue;
    if ( !isShield && fShieldFlag == kShieldOnly ) continue;
    if ( fParticleId != 0 && digihit.pId != fParticleId ) continue;
    
    ndigihit++;
    
    if ( !dshist1 ) dshist1 = hm.Get<TH1F>("dS");
    if ( !dsbyephist1 ) dsbyephist1 = hm.Get<TH1F>("dSByEndPoints");
    if ( !dsbyepvsdsprof) dsbyepvsdsprof=hm.Get<TProfile>("dSByEndPointsvsdS");
    if ( !dehist1 ) dehist1 = hm.Get<TH1F>("dE");
    if ( !dedshist1 ) dedshist1 = hm.Get<TH1F>("dEdS");
    if ( !dedsbyephist1 ) dedsbyephist1 = hm.Get<TH1F>("dEdSByEndPoints");
    if ( !dedsvspeprof) dedsvspeprof = hm.Get<TProfile>("dEdSvspE");
    if ( !dedsbyepvspeprof) dedsbyepvspeprof = hm.Get<TProfile>("dEdSByEndPointsvspE");

    FillTrueDigiHitdS(digihitpath,digihit,dshist1,clause);
    FillTrueDigiHitdSByEndPoints(digihitpath,digihit,dsbyephist1,clause);
    FillTrueDigiHitdSByEndPointsvsdS(digihitpath,digihit,dsbyepvsdsprof,
                                     clause);
    FillTrueDigiHitdE(digihitpath,digihit,dehist1,clause);
    FillTrueDigiHitdEdS(digihitpath,digihit,dedshist1,clause);
    FillTrueDigiHitdEdSByEndPoints(digihitpath,digihit,dedsbyephist1,clause);
    FillTrueDigiHitdEdSvspE(digihitpath,digihit,dedsvspeprof,clause);
    FillTrueDigiHitdEdSByEndPointsvspE(digihitpath,digihit,dedsbyepvspeprof,
                                      clause);

  }

  FillTrueNumDigiHit(digihitpath,ndigihit,clause);
  
  return;
  
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueDigiHitpId(
       std::string digihitpath, const NtpMCDigiScintHit& digihit,
       TH1F* hist1,std::string clause) {
  // Private method to fill ../DigiHit/pId histogram

  if ( !hist1 ) {
    hist1 = new TH1F("pId","DigiHit Particle Id",5000,-2500.,2500.);
    std::string xtitle = "DigiHit Particle Id" + clause;
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/1");
    HistMan hm(digihitpath.c_str());
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(digihit.pId);

}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueDigiHitTrkId(
       std::string digihitpath, const NtpMCDigiScintHit& digihit,
       TH1F* hist1,std::string clause) {
  // Private method to fill ../DigiHit/TrkId histogram

  if ( !hist1 ) {
    hist1 = new TH1F("TrkId","DigiHit StdHep Array Index",200,-100.,100.);
    std::string xtitle = "DigiHit StdHep Array Index" + clause;
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/1");
    HistMan hm(digihitpath.c_str());
    hm.Adopt("",hist1);
  }
  
  hist1 -> Fill(digihit.trkId);

}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueDigiHitdS(
       std::string digihitpath, const NtpMCDigiScintHit& digihit,
       TH1F* hist1,std::string clause) {
  // Private method to fill ../DigiHit/dS histogram

  if ( !hist1 ) {
    hist1 = new TH1F("dS","DigiHits Hit dS",150,0.,5.);
    std::string xtitle = "DigiHit dS (cm)" + clause;
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/0.033 cm");
    HistMan hm(digihitpath.c_str());
    hm.Adopt("",hist1);
  }

  hist1 -> Fill(digihit.dS/Munits::cm);

}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueDigiHitdSByEndPoints(
       std::string digihitpath, const NtpMCDigiScintHit& digihit,
       TH1F* hist1,std::string clause) {
  // Private method to fill ../DigiHit/dS histogram

  if ( !hist1 ) {
    hist1 = new TH1F("dSByEndPoints","DigiHits dS by EndPoints",150,0.,5.);
    std::string xtitle = "DigiHit dS_{linear} (cm)" + clause;
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/0.033 cm");
    HistMan hm(digihitpath.c_str());
    hm.Adopt("",hist1);
  }

  double dx = digihit.x1 - digihit.x0;
  double dy = digihit.y1 - digihit.y0;
  double dz = digihit.z1 - digihit.z0;
  double dsbyendpoints = sqrt(dx*dx+dy*dy+dz*dz)/Munits::cm;

  hist1 -> Fill(dsbyendpoints);

}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueDigiHitdE(
       std::string digihitpath, const NtpMCDigiScintHit& digihit,
       TH1F* hist1,std::string clause) {
  // Private method to fill ../DigiHit/dE

  if ( !hist1 ) {
    hist1 = new TH1F("dE","DigiHits Hit dE",100,0.,5.);
    std::string xtitle = "DigiHit dE (MeV)" + clause;
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/0.05 MeV");
    HistMan hm(digihitpath.c_str());
    hm.Adopt("",hist1);
  }

  hist1 -> Fill(digihit.dE/Munits::MeV);

}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueDigiHitdEdS(
       std::string digihitpath, const NtpMCDigiScintHit& digihit,
       TH1F* hist1,std::string clause) {
  // Private method to fill ../DigiHit/dEdS histogram

  if ( !hist1 ) {
    hist1 = new TH1F("dEdS","DigiHits Hit dE/dS",100,0.,5.);
    std::string xtitle = "DigiHit dE/dS (MeV/cm)" + clause;
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/0.05 MeV/cm");
    HistMan hm(digihitpath.c_str());
    hm.Adopt("",hist1);
  }

  double dEdS = ( digihit.dS > 1.e-10 ) ?
    digihit.dE/Munits::MeV/(digihit.dS/Munits::cm) : 999.;
  
  hist1 -> Fill(dEdS);

}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueDigiHitdEdSByEndPoints(
       std::string digihitpath, const NtpMCDigiScintHit& digihit,
       TH1F* hist1,std::string clause) {
  // Private method to fill ../DigiHit/dEdSByEndPoints histogram

  if ( !hist1 ) {
    hist1 = new TH1F("dEdSByEndPoints","DigiHits Hit dE/dS By EndPoints",
                                                            100,0.,5.);
    std::string xtitle = "DigiHit dE/dS_{linear} (MeV/cm)" + clause;
    hist1->SetXTitle(xtitle.c_str());
    hist1->SetYTitle("Entries/0.05 MeV/cm");
    HistMan hm(digihitpath.c_str());
    hm.Adopt("",hist1);
  }

  double dx = digihit.x1 - digihit.x0;
  double dy = digihit.y1 - digihit.y0;
  double dz = digihit.z1 - digihit.z0;
  double dsbyendpoints = sqrt(dx*dx+dy*dy+dz*dz);

  // for ds greater than a micron
  double dEdS = ( dsbyendpoints > 1.e-10 ) ?
    digihit.dE/Munits::MeV/(dsbyendpoints/Munits::cm) : 999.;
  
  hist1 -> Fill(dEdS);

}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueDigiHitdSByEndPointsvsdS (
       std::string digihitpath, const NtpMCDigiScintHit& digihit,
       TProfile* prof,std::string clause) {
  // Private method to fill ../DigiHit/dEdSvspE profile

  if ( !prof ) { 
    prof = new TProfile("dSByEndPointsvsdS","DigiHit dSByEndPoints vs dS",
                                                    150,0.,5.,0.,5.,"");
    std::string xtitle = "dS (cm)" + clause;
    prof->SetXTitle(xtitle.c_str());
    prof->SetYTitle("dS_{linear} (cm)");
    HistMan hm(digihitpath.c_str());
    hm.Adopt("",prof);
  }

  double dx = digihit.x1 - digihit.x0;
  double dy = digihit.y1 - digihit.y0;
  double dz = digihit.z1 - digihit.z0;
  double dsbyendpoints = sqrt(dx*dx+dy*dy+dz*dz)/Munits::cm;
  prof -> Fill(digihit.dS/Munits::cm,dsbyendpoints);

}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueDigiHitdEdSvspE (
       std::string digihitpath, const NtpMCDigiScintHit& digihit,
       TProfile* prof,std::string clause) {
  // Private method to fill ../DigiHit/dEdSvspE profile

  if ( !prof ) { 
    prof = new TProfile("dEdSvspE","DigiHit dE/dx vs log10(gamma*beta)",
                                                    40,-1.,3.,0.,50.,"");
    std::string xtitle = "log10(#gamma*#beta = p/Mc)" + clause;
    prof->SetXTitle(xtitle.c_str());
    prof->SetYTitle("DigiHit dE/dS (MeV/cm)");
    HistMan hm(digihitpath.c_str());
    hm.Adopt("",prof);
  }

  double dEdS = ( digihit.dS > 1.e-10 ) ?
    (digihit.dE/Munits::MeV)/(digihit.dS/Munits::cm) : -999.;
  if ( dEdS < 0 ) return;
  
  //cout << "dE " << digihit.dE << " dS " << digihit.dS 
  //      << " dEdS (MeV/cm)" << dEdS << endl;
  
  const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
  TParticlePDG* particle = dbpdg.GetParticle(digihit.pId);
  if ( !particle ) {
    MSG("MCMon",Msg::kWarning) << "In FillTrueDigiHitdEdSvspE:\n"
                               << "No particle of type " << digihit.pId
                               << " in PDG database. Skipped hit." << endl;
    return;
  }
  
  Double_t partmass = particle->Mass();
  //cout << "partmass " << partmass << endl;
  //cout << "digihit.pE " << digihit.pE << endl;
 
  Double_t betagamma2 = pow(digihit.pE/partmass,2) - 1.;
  if ( betagamma2 < 0. ) return; // precision errors
  
  Double_t log10betagamma = log10(sqrt(betagamma2));
  //cout << "log10(beta*gamma) " << log10betagamma << endl;
  
  prof -> Fill(log10betagamma,dEdS);

}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::FillTrueDigiHitdEdSByEndPointsvspE (
       std::string digihitpath, const NtpMCDigiScintHit& digihit,
       TProfile* prof,std::string clause) {
  // Private method to fill ../DigiHit/dEdSByEndPointsvspE profile

  if ( !prof ) { 
    prof = new TProfile("dEdSByEndPointsvspE",
                        "DigiHit dE/dx_{linear} vs log10(gamma*beta)",
                                                    40,-1.,3.,0.,50.,"");
    std::string xtitle = "log10(#gamma*#beta = p/Mc)" + clause;
    prof->SetXTitle(xtitle.c_str());
    prof->SetYTitle("DigiHit dE/dS (MeV/cm)");
    HistMan hm(digihitpath.c_str());
    hm.Adopt("",prof);
  }

  double dx = digihit.x1 - digihit.x0;
  double dy = digihit.y1 - digihit.y0;
  double dz = digihit.z1 - digihit.z0;
  double dsbyendpoints = sqrt(dx*dx+dy*dy+dz*dz);

  double dEdS = ( dsbyendpoints > 1.e-10 ) ?
    (digihit.dE/Munits::MeV)/(dsbyendpoints/Munits::cm) : -999.;
  if ( dEdS < 0 ) return;
  
  
  const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
  TParticlePDG* particle = dbpdg.GetParticle(digihit.pId);
  if ( !particle ) {
    MSG("MCMon",Msg::kWarning) << "In FillTrueDigiHitdEdSByEndPointsvspE:\n"
                               << "No particle of type " << digihit.pId
                               << " in PDG database. Skipped hit." << endl;
    return;
  }
  
  Double_t partmass = particle->Mass();
 
  Double_t betagamma2 = pow(digihit.pE/partmass,2) - 1.;
  if ( betagamma2 < 0. ) return; // precision errors
  
  Double_t log10betagamma = log10(sqrt(betagamma2));
  
  prof -> Fill(log10betagamma,dEdS);

}

  
//_______________________________________________________________________
int MCMonitorCosmicHistograms::NumHistograms(const TFolder* baseFolder) const {
  // Default argument null => start at base folder
  // number of histograms in the class
  
  int numHist = 0;
  
  if ( !baseFolder ) {
    HistMan hm(fName.c_str());
    const TFolder* folder = dynamic_cast<TFolder*>
                            (hm.BaseFolder().FindObject(fName.c_str()));
    if ( folder ) return NumHistograms(folder);
  }
  else {
    cout << "In NumHistograms for folder " << baseFolder -> GetName() << endl;

    TCollection* objlist = baseFolder -> GetListOfFolders();
    TIter nextobj(objlist);
    TObject* obj = 0;
    while ( ( obj = (TObject*)nextobj() ) ) {
      cout << "Printing object" << endl;
      obj -> Print();
      TFolder* folder = dynamic_cast<TFolder*>(obj);
      if ( folder ) numHist += this -> NumHistograms(folder);
      else numHist++;
    }
  }
  
  return numHist;
   
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::SetCanRebin(const TFolder* baseFolder) {
  // Default argument null => start at base folder
  // Set all the histograms so they can be rebinned by the user by default

  if ( !baseFolder ) {
    HistMan hm(fName.c_str());
    const TFolder* folder = dynamic_cast<TFolder*>
                            (hm.BaseFolder().FindObject(fName.c_str()));
    if ( folder ) return SetCanRebin(folder);
  }
  else {
    cout << "In SetCanRebin for folder " << baseFolder -> GetName() << endl;

    TCollection* objlist = baseFolder -> GetListOfFolders();
    TIter nextobj(objlist);
    TObject* obj = 0;
    while ( ( obj = (TObject*)nextobj() ) ) {
      cout << "Printing object" << endl;
      obj -> Print();
      TFolder* folder = dynamic_cast<TFolder*>(obj);
      if ( folder ) this -> SetCanRebin(folder);
      else {
        TH1* hist = dynamic_cast<TH1*>(obj);
        hist -> SetBit(TH1::kCanRebin);
      }
    }
  }
  
  return;
   
}

//_______________________________________________________________________
void MCMonitorCosmicHistograms::AllSumw2(const TFolder* baseFolder) {
  // Default argument null => start at base folder
  // sets default Poisson errors on all histograms
  
  if ( !baseFolder ) {
    HistMan hm(fName.c_str());
    const TFolder* folder = dynamic_cast<TFolder*>
                            (hm.BaseFolder().FindObject(fName.c_str()));
    if ( folder ) return AllSumw2(folder);
  }
  else {
    cout << "In AllSumw2 for folder " << baseFolder -> GetName() << endl;

    TCollection* objlist = baseFolder -> GetListOfFolders();
    TIter nextobj(objlist);
    TObject* obj = 0;
    while ( ( obj = (TObject*)nextobj() ) ) {
      cout << "Printing object" << endl;
      obj -> Print();
      TFolder* folder = dynamic_cast<TFolder*>(obj);
      if ( folder ) this -> AllSumw2(folder);
      else {
        TH1* hist = dynamic_cast<TH1*>(obj);
        if ( hist ) {
          TProfile* prof = dynamic_cast<TProfile*>(obj);
          if (!prof) hist -> Sumw2();
        }
      }
    }
  }
  
  return;
   
}

---