loon: ANtpInfoObjectFillerMRCC.cxx Source File

00001 #include "ANtpInfoObjectFillerMRCC.h"
00002 
00003 #include "../ANtpEventInfoMRCC.h"
00004 #include "../ANtpTruthInfoBeam.h"
00005 
00006 #include "ANtpInfoObjectFiller.h"
00007 #include "ANtpInfoObjectFillerBeam.h"
00008 
00009 #include "MuonRemoval/NtpMREvent.h"
00010 #include "MuonRemoval/NtpMRRecord.h"
00011 #include "StandardNtuple/NtpStRecord.h"
00012 #include "Conventions/DetectorType.h"
00013 #include "Conventions/BeamType.h"
00014 #include "MessageService/MsgService.h"
00015 #include "Mad/MadAbID.h"
00016 #include "PhysicsNtuple/Store/Interface.h"
00017 
00018 #include "TClonesArray.h"
00019 #include "TSystem.h"
00020 
00021 #include <cassert>
00022 #include <string>
00023 #include <vector>
00024 
00025 using namespace std;
00026 
00027 CVSID("$Id: ANtpInfoObjectFillerMRCC.cxx,v 1.3 2009/11/24 16:50:39 rodriges Exp $");
00028 
00029 void ANtpInfoObjectFillerMRCC::
00030 FillEventInformation(int i,
00031                      const NtpMRRecord* mrRecord,
00032                      const NtpStRecord* stRecord,
00033                      ANtpEventInfoMRCC* antp,
00034                      ANtpTruthInfoBeam* truth) const
00035 {
00036   assert(i >= 0);
00037   assert(mrRecord);
00038 
00039   // The number of mrcc records in this snarl
00040   const int J = mrRecord->mrevt->GetEntriesFast();
00041   // Go through all the mrcc records in this snarl
00042   for(int j = 0; j < J; ++j){
00043     const NtpMREvent* mrevt = (NtpMREvent*)(mrRecord->mrevt->At(j));
00044     assert(mrevt);
00045     // See if any of them are matched up to the current event
00046     if(mrevt->best_event == i){
00047       FillEventInformationHelper(mrevt, stRecord, antp, truth);
00048       // Assume there is only one match. If there are more then there's
00049       // nothing sensible we can do about it anyway.
00050       return;
00051     }
00052   } // end for j
00053   // No matches, leave antp alone (it should already be default)
00054 }
00055 
00056 
00057 void ANtpInfoObjectFillerMRCC::
00058 FillEventInformationHelper(const NtpMREvent* mrnt,
00059                            const NtpStRecord* sntp_orig,
00060                            ANtpEventInfoMRCC* antp,
00061                            ANtpTruthInfoBeam* truth) const
00062 {
00063   assert(mrnt);
00064   assert(sntp_orig);
00065   assert(antp);
00066 
00067   antp->index = mrnt->index;
00068   antp->digits = mrnt->ndigit;
00069   antp->strips = mrnt->nstrip;
00070 
00071   antp->origEvent = mrnt->orig_event;
00072   antp->bestEvent = mrnt->best_event;
00073   antp->bestPurity = mrnt->best_purity;
00074   antp->bestComplete = mrnt->best_complete;
00075   antp->bestPurityPHw = mrnt->best_purity_phw;
00076   antp->bestCompletePHw = mrnt->best_complete_phw;
00077 
00078   MSG("ObjectFillerMRCC",Msg::kDebug) 
00079     << "orig_event=" << mrnt->orig_event << endl
00080     //<< "mrnt->snarl=" << mrnt->GetHeader().GetSnarl() << endl
00081     << "sntp_orig->snarl=" << sntp_orig->GetHeader().GetSnarl() << endl;
00082 
00083   vector<const NtpSREvent*> orig_evts=sntp_orig->GetEvents();
00084   assert(mrnt->orig_event < (int)orig_evts.size());
00085 
00086   const NtpSREvent* orig_event=orig_evts[mrnt->orig_event];
00087   antp->origABPID = ABPID(orig_event, sntp_orig);
00088   antp->origROPID = ROPID(orig_event, sntp_orig);
00089 
00090   if(truth){
00091     // Fill the truth information. This code copied from
00092     // ANtpInfoObjectFillerNC::FillMCInformation
00093     truth->Reset();
00094     // Gah, const correctness, or lack of it
00095     NtpMCTruth* ntpMCTruth=const_cast<NtpMCTruth*>(sntp_orig->GetEventMCTruth(mrnt->orig_event));
00096     
00097     ANtpInfoObjectFiller::FillMCTruthInformation(ntpMCTruth, truth);
00098     ANtpInfoObjectFillerBeam::FillBeamMCTruthInformation(ntpMCTruth, sntp_orig->stdhep, 
00099                                                          truth);
00100   }
00101 
00102   antp->showerVtxX = mrnt->shwvtxx;
00103   antp->showerVtxY = mrnt->shwvtxy;
00104   antp->showerVtxZ = mrnt->shwvtxz;
00105   antp->showerEndX = mrnt->shwendx;
00106   antp->showerEndY = mrnt->shwendy;
00107   antp->showerEndZ = mrnt->shwendz;
00108   antp->showerVtxPlane = mrnt->shwvtxplane;
00109   antp->showerEndPlane = mrnt->shwendplane;
00110   antp->showerPlanes = mrnt->shwnplane;
00111   antp->showerCharge = mrnt->shwcharge;
00112 
00113   antp->trackVtxX = mrnt->vtxx;
00114   antp->trackVtxY = mrnt->vtxy;
00115   antp->trackVtxZ = mrnt->vtxz;
00116   antp->trackVtxDistanceToEdge = mrnt->vtxdistance;
00117   antp->trackEndX = mrnt->endx;
00118   antp->trackEndY = mrnt->endy;
00119   antp->trackEndZ = mrnt->endz;
00120   antp->trackEndDistanceToEdge = mrnt->enddistance;
00121 
00122   antp->trackVtxPlane = mrnt->vtxp;
00123   antp->trackEndPlane = mrnt->endp;
00124   antp->trackPlanes = mrnt->npln;
00125   antp->rangeMomentum = mrnt->prng;
00126   antp->fitMomentum = mrnt->pcrv;
00127   antp->trackDcosXVtx = mrnt->pvdx;
00128   antp->trackDcosYVtx = mrnt->pvdy;
00129   antp->trackDcosZVtx = mrnt->pvdz;
00130 
00131   antp->passedFit = mrnt->fitp;
00132   antp->contained = mrnt->endc;
00133   antp->containedOrGoodFit = mrnt->pass;
00134   antp->trackMomentumX = mrnt->pmux;
00135   antp->trackMomentumY = mrnt->pmuy;
00136   antp->trackMomentumZ = mrnt->pmuz;
00137   antp->maxTrackPlane = mrnt->mxpl;
00138   antp->qp = mrnt->qp;
00139 
00140   antp->trueMomentumX = mrnt->mrmpmux;
00141   antp->trueMomentumY = mrnt->mrmpmuy;
00142   antp->trueMomentumZ = mrnt->mrmpmuz;
00143   antp->trueQ2 = mrnt->mrmQ2;
00144   antp->trueShowerEnergy = mrnt->mrmEshw;
00145 
00146 }
00147 
00148 double ANtpInfoObjectFillerMRCC::
00149 ABPID(const NtpSREvent* ev, const NtpStRecord* st) const
00150 {
00151   static bool once=false;
00152   static MadAbID ABID;
00153   
00154   // Code stolen from NueAna's MuonRemovalInfoAna
00155   if(!once){
00156     once=true;
00157     string base=getenv("SRT_PRIVATE_CONTEXT");
00158      if(base!="" && base!="."){
00159        // check if directory exists in SRT_PRIVATE_CONTEXT
00160        string path = base + "/Mad/data";
00161        void* dir_ptr = gSystem->OpenDirectory(path.c_str());
00162        if(!dir_ptr){
00163          base=getenv("SRT_PUBLIC_CONTEXT");
00164        }
00165      }
00166      else{  
00167        base=getenv("SRT_PUBLIC_CONTEXT"); 
00168      }
00169      if(base=="") {
00170        MSG("ObjectFillerMRCC",Msg::kFatal)<<"No SRT_PUBLIC_CONTEXT set "
00171                                          <<"Do not know where to look "
00172                                          <<"for AB pdf files "<<std::endl;
00173        // Dummy message to make sure it exits
00174        MSG("ObjectFillerMRCC",Msg::kFatal) << "Goodbye" << endl;
00175      }
00176      base+="/Mad/data";
00177      string fname=base;
00178      string rmc="";
00179      // ARGH! Can't do this at all sensibly, since there aren't PDF
00180      // files for dogwood (I guess). Should also really check on MC
00181      // version, but this code is never going to be run on anything
00182      // other than daikon, so it doesn't matter.
00183      ReleaseType::Release_t release=st->GetRelease();
00184      if(ReleaseType::IsCedar(release)){
00185        rmc="cedar_daikon";
00186      }
00187      else{
00188        MSG("ObjectFillerMRCC",Msg::kWarning) 
00189          << "Running on release type " << ReleaseType::AsString(release)
00190          << " which has no PDFs. Using cedar_daikon instead." << endl;
00191        rmc="cedar_daikon";
00192      }
00193 
00194      string sdet="";
00195      Detector::Detector_t det=st->GetVldContext()->GetDetector();
00196      if(det==Detector::kNear){
00197        sdet="near";
00198      }
00199      else if(det==Detector::kFar){
00200        sdet="far";
00201      }
00202      else{
00203        MSG("ObjectFillerMRCC",Msg::kFatal) << "Dont know detector type."
00204                                            << " Bailing"<<endl;
00205        // Dummy message to make sure it exits
00206        MSG("ObjectFillerMRCC",Msg::kFatal) << "Goodbye" << endl;
00207      }
00208      string sbeam="le";
00209      /*
00210        Should really do this switch, but I don't know how to get the
00211        beamtype through to here to check against. Probably won't ever
00212        run it on non-LE events anyway.
00213 
00214      switch(current_beam){
00215      case BeamType::kL010z000i:
00216        sbeam="le0";
00217        break;
00218      case BeamType::kL010z170i:
00219        sbeam="le170";
00220        break;
00221      case BeamType::kL010z185i:
00222        sbeam="le";
00223        break;
00224      case BeamType::kL010z200i:
00225        sbeam="le200";
00226        break;
00227      case BeamType::kL100z200i:
00228        sbeam="pme";
00229        break;
00230      case BeamType::kL150z200i:
00231        sbeam="pme";
00232        break;
00233      case BeamType::kL250z200i:
00234        sbeam="phe";
00235        break;
00236      case BeamType::kLE:
00237        sbeam="le";
00238        break;
00239      default:
00240        MSG("ObjectFillerMRCC",Msg::kWarning)<<"Don't know beam type "
00241                                            <<" defaulting to LE"<<endl;
00242        sbeam="le";
00243        break;
00244      }
00245      */
00246                                                                                  
00247      fname+="/ab_pdf_"+sdet+"_"+sbeam+"_"+rmc+".root";
00248                                                                                  
00249      ABID.ReadPDFs(fname.c_str());
00250   }
00251                                                                                  
00252   return ABID.CalcPID(ev, st);
00253 }
00254 
00255 double ANtpInfoObjectFillerMRCC::
00256 ROPID(const NtpSREvent* /*ev*/, const NtpStRecord* /*st*/) const
00257 {
00258   // ***********
00259   // Don't calculate the kNN stuff. It's too slow
00260   // ***********
00261   return -1;
00262 
00263   /*
00264   // TODO: This code copied wholesale from ANtpInfoObjectFillerNC. Should rationalize.
00265 
00266   static bool once=false;
00267   static Anp::Interface anpint;
00268 
00269   if(!once){
00270     once=true;
00271     //set the location of the weight files based on the detector
00272     TString weightFileName = "knn.train.far.cedar.daikon.root";
00273     Detector::Detector_t det=st->GetVldContext()->GetDetector();
00274     if(det == Detector::kNear) weightFileName = "knn.train.near.cedar.daikon.root";
00275 
00276     TString base = getenv("SRT_PRIVATE_CONTEXT");
00277     TString ncutils = "/NCUtils/data/";
00278     if(base != "" && base != "."){
00279       //check that the private context director exists, if not use public context
00280       void *dirptr = gSystem->OpenDirectory(base+ncutils);
00281       if(!dirptr){
00282         base = getenv("SRT_PUBLIC_CONTEXT");
00283       }//end if private context doesnt exist
00284     }//end if private context is defined
00285     else base = getenv("SRT_PUBLIC_CONTEXT");
00286 
00287     //check that the public context exists, if not assert false
00288     if(base == ""){
00289       MSG("ObjectFillerMRCC", Msg::kFatal) << "no SRT_PUBLIC_CONTEXT set"
00290                                                  << endl;
00291       assert(false);
00292     }
00293 
00294     TString weightFilePath = base+ncutils+weightFileName;
00295     TString baseConf = getenv("SRT_PUBLIC_CONTEXT");
00296 
00297     Registry reg;
00298     reg.Set("InterfaceConfigPath", 
00299             baseConf+"/PhysicsNtuple/Config/Config2007Real.txt");
00300     reg.Set("FillkNNFilePath", weightFilePath);
00301     
00302     anpint.Config(reg);
00303   }
00304 
00305   static int prevSnarl=-1;
00306   int thisSnarl=st->GetHeader().GetSnarl();
00307   // Only refill the Anp interface once per snarl
00308   if(thisSnarl!=prevSnarl){
00309     prevSnarl=thisSnarl;
00310     // Gah, I wish const correctness wasn't such a pain
00311     if(!anpint.FillSnarl(const_cast<NtpStRecord*>(st))){
00312       MSG("ObjectFillerMRCC", Msg::kFatal) << "Couldn't set up Anp::Interface" << endl;
00313       // Dummy message to make sure it exits here
00314       MSG("ObjectFillerMRCC", Msg::kFatal) << "Goodbye" << endl;
00315     }
00316   }
00317 
00318   // Gah, I wish const correctness wasn't such a pain
00319   return anpint.GetVar("knn_pid", const_cast<NtpSREvent*>(ev));
00320   */
00321 }