loon: NuPIDInterface.cxx Source File

00001 
00004 // Coded by Jeff Hartnell Jan/2007 onwards
00005 //
00006 // Contact: j.j.hartnell@rl.ac.uk
00008 
00009 #include <fstream>
00010 
00011 #include "Conventions/Detector.h"
00012 // NuBarPID Officially Depracated
00013 //#include "NuBarPID/LoadPDF.h"
00014 #include "Mad/MadAbID.h"
00015 #include "Mad/MadDpID.h"
00016 #include "MessageService/MsgService.h"
00017 #include "NtupleUtils/NuEvent.h"
00018 #include "NtupleUtils/NuLibrary.h"
00019 // NuBarPID Officially Depracated
00020 //#include "NuBarPID/PIDCalculator.h"
00021 #include "Registry/Registry.h"
00022 #include "Conventions/ReleaseType.h"//used here from 25/Feb/07
00023 
00024 #include "NuMuBar/NuMadAnalysis.h"
00025 #include "NuMuBar/NuPIDInterface.h"
00026 
00027 using std::endl;
00028 using std::cout;
00029 
00030 CVSID("$Id: NuPIDInterface.cxx,v 1.31 2009/12/08 18:13:00 ahimmel Exp $");
00031 
00032 //......................................................................
00033 
00034 NuPIDInterface::NuPIDInterface()
00035 {
00036   MSG("NuPIDInterface",Msg::kDebug)
00037     <<"Running NuPIDInterface Constructor..."<<endl;
00038 
00039   //initialise pointers
00040   poLoadPDF=0;
00041   poLoadPDFKin=0;
00042   poPIDCalculator=0;
00043   poPIDCalculatorKin=0;
00044   abID=0;
00045   dpID=0;
00046   mad=0;
00047   roInterface=0;
00048   jmInterface=0;
00049   roRegistry=0;
00050   
00051   MSG("NuPIDInterface",Msg::kDebug)
00052     <<"Finished NuPIDInterface Constructor"<<endl;
00053 }
00054 
00055 //......................................................................
00056 
00057 NuPIDInterface::~NuPIDInterface()
00058 {
00059   MSG("NuPIDInterface",Msg::kDebug)
00060     <<"Running NuPIDInterface Destructor..."<<endl;
00061 
00062   //if (poLoadPDF) delete poLoadPDF;
00063   //if (poLoadPDFKin) delete poLoadPDFKin;
00064   //if (poPIDCalculator) delete poPIDCalculator;
00065   //if (poPIDCalculatorKin) delete poPIDCalculatorKin;
00066   //if (abID) delete abID;
00067   //if (dpID) delete dpID;
00068   //if (mad) delete mad;
00069   //if (roInterface) delete roInterface;
00070   //if (roRegistry) delete roRegistry;
00071 
00072   MSG("NuPIDInterface",Msg::kDebug)
00073     <<"Finished NuPIDInterface Destructor"<<endl;
00074 }
00075 
00076 //......................................................................
00077 // PoID Officially depracated
00078 //std::string NuPIDInterface::GetFileNamePoID(const NuEvent& nu) const
00079 //{
00080 //  //get an instance of the code library
00081 //  const NuLibrary& lib=NuLibrary::Instance();
00082 //
00083 //  string sBase=lib.general.GetReleaseDirToUse("NuBarPID");
00084 //  string sFileName=sBase+"/NuBarPID/files/POpdf_birch_carrot.root";
00085 //  
00086 //  Bool_t badVersions=false;
00087 //  
00088 //  //there is just one file to choose from as of 20070816
00089 //  sFileName=sBase+"/NuBarPID/files/POpdf_cedar_daikon_ND.root";
00090 //
00092 //  if (nu.detector==Detector::kNear) {
00093 //    if (ReleaseType::IsBirch(nu.recoVersion) &&
00094 //      ReleaseType::IsCarrot(nu.mcVersion)) {
00095 //      sFileName=sBase+"/NuBarPID/files/POpdf_birch_carrot.root";
00096 //    }
00097 //    else if (ReleaseType::IsBirch(nu.recoVersion) &&
00098 //           ReleaseType::IsDaikon(nu.mcVersion)) {
00099 //      sFileName=sBase+"/NuBarPID/files/POpdf_birch_pre-daikon.root";
00100 //      MSG("NuPIDInterface",Msg::kWarning)
00101 //      <<"For Birch-Daikon using the wrong filename="<<sFileName<<endl;
00102 //    }
00103 //    else if (ReleaseType::IsCedar(nu.recoVersion) &&
00104 //           ReleaseType::IsCarrot(nu.mcVersion)) {
00105 //      sFileName=sBase+"/NuBarPID/files/POpdf_cedar_carrot.root";
00106 //    }
00107 //    else if (ReleaseType::IsCedar(nu.recoVersion) &&
00108 //           ReleaseType::IsDaikon(nu.mcVersion)) {
00109 //      sFileName=sBase+"/NuBarPID/files/POpdf_cedar_daikon_ND.root";
00110 //    }
00111 //    else {
00112 //      badVersions=true;
00113 //    }
00114 //  }
00115 //  else if (nu.detector==Detector::kFar) {
00116 //    if (ReleaseType::IsBirch(nu.recoVersion) &&
00117 //      ReleaseType::IsCarrot(nu.mcVersion)) {
00118 //      sFileName=sBase+"/NuBarPID/files/POpdf_birch_carrot_FD.root";
00119 //    }
00120 //    else if (ReleaseType::IsBirch(nu.recoVersion) &&
00121 //           ReleaseType::IsDaikon(nu.mcVersion)) {
00122 //      sFileName=sBase+"/NuBarPID/files/POpdf_birch_carrot_FD.root";
00123 //      MSG("NuPIDInterface",Msg::kWarning)
00124 //      <<"For Birch-Daikon using the wrong filename="<<sFileName<<endl;
00125 //    }
00126 //    else if (ReleaseType::IsCedar(nu.recoVersion) &&
00127 //           ReleaseType::IsCarrot(nu.mcVersion)) {
00128 //      sFileName=sBase+"/NuBarPID/files/POpdf_cedar_carrot_FD.root";
00129 //    }
00130 //    else if (ReleaseType::IsCedar(nu.recoVersion) &&
00131 //           ReleaseType::IsDaikon(nu.mcVersion)) {
00132 //      sFileName=sBase+"/NuBarPID/files/POpdf_cedar_carrot_FD.root";
00133 //      MSG("NuPIDInterface",Msg::kWarning)
00134 //      <<"For Cedar-Daikon using the wrong filename="<<sFileName<<endl;
00135 //    }
00136 //    else {
00137 //      badVersions=true;
00138 //    }
00139 //  }
00140 //  else {
00141 //    cout<<"Ahhh, detector="<<nu.detector<<" not recognised"<<endl;
00142 //  }
00143 //*/
00144 //
00145 //  if (badVersions) {
00146 //    MSG("NuPIDInterface",Msg::kError)
00147 //      <<"Ahhh, ReleaseType not recognised:"<<endl;
00148 //    this->PrintReleaseTypeEtc(nu);
00149 //    MSG("NuPIDInterface",Msg::kError)<<"Will abort here..."<<endl;
00150 //    assert(false);
00151 //  }
00152 //  
00153 //  return sFileName;
00154 //}
00155 
00156 //......................................................................
00157 // PoID Officially depracated
00158 //std::string NuPIDInterface::GetFileNamePoIDKin(const NuEvent& nu) const
00159 //{
00160 //  //get an instance of the code library
00161 //  const NuLibrary& lib=NuLibrary::Instance();
00162 //
00163 //  string sBase=lib.general.GetReleaseDirToUse("NuBarPID");
00164 //  string sFileName=sBase+"/NuBarPID/files/POpdf_birch_carrot.root";
00165 //  
00166 //  Bool_t badVersions=false;
00167 //  
00168 //  //there is just one file to choose from as of 20080222
00169 //  sFileName=sBase+"/NuBarPID/files/POpdf_cedar_daikon_ND_Kin.root";
00170 //
00171 //  if (badVersions) {
00172 //    MSG("NuPIDInterface",Msg::kError)
00173 //      <<"Ahhh, ReleaseType not recognised:"<<endl;
00174 //    this->PrintReleaseTypeEtc(nu);
00175 //    MSG("NuPIDInterface",Msg::kError)<<"Will abort here..."<<endl;
00176 //    assert(false);
00177 //  }
00178 //  
00179 //  return sFileName;
00180 //}
00181 
00182 //......................................................................
00183 
00184 std::string NuPIDInterface::GetFileNameAbID(const NuEvent& nu) const
00185 {
00186   //get an instance of the code library
00187   const NuLibrary& lib=NuLibrary::Instance();
00188   
00189   string sBase=lib.general.GetReleaseDirToUse("Mad");
00190   string sFileName=sBase+"/Mad/data/ab_pdf_near_le_cedar_daikon.root";
00191   Bool_t badVersions=false;
00192 
00193   //Set it to defaults that don't exist so the code stops if we don't
00194   //set the release type
00195   string sVegetable = "Aubergine";
00196   string sTree = "Bonsai";
00197   string sBeam = "LHC";
00198   string sDetector = "Nova";
00199 
00200   // Assuming cedar since these are the only available training files
00201   sTree = "cedar";
00202   if (!ReleaseType::IsCedar(nu.recoVersion)){
00203     MAXMSG("NuPIDInterface",Msg::kInfo,1) << "Assuming cedar reco version in AbID initialization (changed from "
00204     << NuUtilities::PrintRelease(nu.releaseType) << ") since those are the only available training files" << endl;
00205   }
00206 
00207   if (ReleaseType::IsDaikon(nu.mcVersion)){
00208     sVegetable = "daikon";
00209   }
00210   else{
00211     badVersions=true;
00212   }
00213   
00214   if (nu.detector == Detector::kNear){
00215     sDetector = "near";
00216     if (nu.beamType==BeamType::kL010z185i ||
00217         nu.beamType==BeamType::kL010z185i_rev) {
00218       sBeam = "le";
00219     }
00220     else if (nu.beamType==BeamType::kL010z170i) {
00221       sBeam = "le170";
00222     }
00223     else if (nu.beamType==BeamType::kL010z200i) {
00224       sBeam = "le200";
00225     }
00226     else if (nu.beamType==BeamType::kL010z000i) {
00227       sBeam = "le0";
00228     }
00229     else if (nu.beamType==BeamType::kL100z200i ||
00230              nu.beamType==BeamType::kL150z200i ||
00231              nu.beamType==BeamType::kL200z200i) {
00232       //is this correct? Not sure what pme really is
00233       sBeam = "pme";
00234     }
00235     else if (nu.beamType==BeamType::kL250z200i) {
00236       sBeam = "phe";
00237     }
00238     else {
00239       MSG("NuPIDInterface",Msg::kError)
00240         <<"Ahhh, BeamType not recognised: "<<nu.beamType<<endl;
00241       assert(false);
00242     }
00243   }
00244   else if (nu.detector==Detector::kFar) {
00245     sDetector = "far";
00246     //check if beamtype is phe
00247     if (nu.beamType==BeamType::kL250z200i) {
00248       sBeam = "phe";
00249     }
00250     else {
00251       sBeam = "le";
00252     }
00253   }
00254   else {
00255     cout<<"Ahhh, detector="<<nu.detector<<" not recognised"<<endl;
00256   }
00257   
00258   sFileName = sBase
00259     + "/Mad/data/ab_pdf_"
00260     + sDetector + "_"
00261     + sBeam + "_"
00262     + sTree + "_"
00263     + sVegetable + ".root";
00264   /*
00265   if (nu.detector==Detector::kNear) {
00266     if (ReleaseType::IsCedar(nu.recoVersion) &&
00267         ReleaseType::IsDaikon(nu.mcVersion)) {
00268       //check the beam type
00269       if (nu.beamType==BeamType::kL010z185i) {
00270         sFileName=sBase+"/Mad/data/ab_pdf_near_le_cedar_daikon.root";
00271       }
00272       else if (nu.beamType==BeamType::kL010z170i) {
00273         sFileName=sBase+
00274           "/Mad/data/ab_pdf_near_le170_cedar_daikon.root";
00275       }
00276       else if (nu.beamType==BeamType::kL010z200i) {
00277         sFileName=sBase+
00278           "/Mad/data/ab_pdf_near_le200_cedar_daikon.root";
00279       }
00280       else if (nu.beamType==BeamType::kL010z000i) {
00281         sFileName=sBase+
00282           "/Mad/data/ab_pdf_near_le0_cedar_daikon.root";
00283       }
00284       else if (nu.beamType==BeamType::kL100z200i ||
00285                nu.beamType==BeamType::kL150z200i ||
00286                nu.beamType==BeamType::kL200z200i) {
00287         //is this correct? Not sure what pme really is
00288         sFileName=sBase+
00289           "/Mad/data/ab_pdf_near_pme_cedar_daikon.root";
00290       }
00291       else if (nu.beamType==BeamType::kL250z200i) {
00292         sFileName=sBase+
00293           "/Mad/data/ab_pdf_near_phe_cedar_daikon.root";
00294       }
00295       else {
00296         MSG("NuPIDInterface",Msg::kError)
00297           <<"Ahhh, BeamType not recognised: "<<nu.beamType<<endl;
00298         assert(false);
00299       }
00300     }
00301     else {
00302       badVersions=true;
00303     }
00304   }
00305   else if (nu.detector==Detector::kFar) {
00306     if (ReleaseType::IsCedar(nu.recoVersion) &&
00307         ReleaseType::IsDaikon(nu.mcVersion)) {
00308       //check if beamtype is phe
00309       if (nu.beamType==BeamType::kL250z200i) {
00310         sFileName=sBase+"/Mad/data/ab_pdf_far_phe_cedar_daikon.root";
00311       }
00312       else {
00313         sFileName=sBase+"/Mad/data/ab_pdf_far_le_cedar_daikon.root";    
00314       }
00315     }
00316     else {
00317       badVersions=true;
00318     }
00319   }
00320   else {
00321     cout<<"Ahhh, detector="<<nu.detector<<" not recognised"<<endl;
00322   }
00323   */
00324   if (badVersions) {
00325     MSG("NuPIDInterface",Msg::kError)
00326       <<"Ahhh, ReleaseType not recognised:"<<endl;
00327     this->PrintReleaseTypeEtc(nu);
00328     MSG("NuPIDInterface",Msg::kError)<<"Will abort here..."<<endl;
00329     assert(false);
00330   }
00331 
00332   return sFileName;
00333 }
00334 
00335 //......................................................................
00336 std::string NuPIDInterface::GetFileNameRoID(const NuEvent& nu) const
00337 {
00338   return GetFileNamekNNID(nu,1);
00339 }
00340 std::string NuPIDInterface::GetFileNameJmID(const NuEvent& nu) const 
00341 {
00342   return GetFileNamekNNID(nu,3);
00343 }
00344 std::string NuPIDInterface::GetFileNamekNNID(const NuEvent& nu, int type=1) const
00345 {
00346 
00347   //get an instance of the code library
00348   const NuLibrary& lib=NuLibrary::Instance();
00349 
00351   //NOTE: the input files are no longer distributed so they 
00352   //have to be copied to the local site
00353   //the NuMuBar/data directory is the place to copy or soft
00354   //link the files to
00356   
00357   string sBase=lib.general.GetReleaseDirToUse("NuMuBar");
00358   //string sFileName=sBase+"/Mad/data/roid.pid.near.cedar.daikon.root";
00359   string sFileName=sBase+"/NuMuBar/data/knn.physics.near.daikon_04.dogwood1.L010z185i.root";
00360   if(type ==2) sFileName=sBase+"/NuMuBar/data/knn.physics.20p3GeV.near.daikon_04.dogwood1.L010z185i.root";
00361   if(type ==3) sFileName=sBase+"/NuMuBar/data/knn.physics.Alt.near.daikon_04.dogwood1.L010z185i.root";
00362 
00363   Bool_t badVersions=false;
00364 
00365   //Set it to defaults that don't exist so the code stops if we don't
00366   //set the release type
00367   string sVegetable = "Aubergine";
00368   string sTree = "Bonsai";
00369   string sBeam = "LHC";
00370   string sDetector = "Nova";
00371 
00372 //  if (ReleaseType::IsCedar(nu.recoVersion)){
00373 //    //
00374 //    // Cedar files for "Alt" do not exist 
00375 //    //
00376 //    sTree = "dogwood1";
00377 //  }
00378 //  else if (ReleaseType::IsDogwood(nu.recoVersion)){
00379 //    sTree = "dogwood1";
00380 //  }
00381 //  else {
00382 //    badVersions=true;
00383 //  }
00384   // Assuming dogwood1 since these are the only available training files
00385   sTree = "dogwood1";
00386   if (!ReleaseType::IsDogwood(nu.recoVersion)){
00387     MAXMSG("NuPIDInterface",Msg::kInfo,1) << "Assuming dogwood1 reco version in RoID initialization (changed from "
00388     << NuUtilities::PrintRelease(nu.releaseType) << ") since those are the only available training files" << endl;
00389   }
00390   
00391 
00392   if (ReleaseType::IsDaikon(nu.mcVersion)){
00393     sVegetable = "daikon_04";
00394   }
00395   else{
00396     badVersions=true;
00397   }
00398 
00399   if (nu.detector==Detector::kNear) {
00400     sDetector = "near";
00401   }
00402   else if (nu.detector==Detector::kFar) {
00403     sDetector = "far";
00404   }
00405   else {
00406     cout<<"Ahhh, detector="<<nu.detector<<" not recognised"<<endl;
00407   }
00408 
00409   //There is no alternative:
00410   sBeam = "L010z185i";
00411   if(type==1)
00412     {
00413       sFileName = sBase + "/NuMuBar/data/knn.physics.";
00414     }
00415   else if(type==2)
00416     { 
00417       sFileName = sBase + "/NuMuBar/data/knn.physics.20p3GeV.";
00418     }
00419   else if(type==3)
00420     {
00421       sFileName = sBase + "/NuMuBar/data/knn.physics.Alt.";
00422     }
00423   else{
00424     cerr<<"NuPIDInterface::GetFileNamekNNID  -- invalid type "<<type<<endl;
00425     assert(false);
00426   }
00427   sFileName+= sDetector + "."
00428     + sVegetable + "."
00429     + sTree + "."
00430     + sBeam + ".root";
00431 
00432   if (badVersions) {
00433     MSG("NuPIDInterface",Msg::kError)
00434       <<"Ahhh, ReleaseType not recognised:"<<endl;
00435     this->PrintReleaseTypeEtc(nu);
00436     MSG("NuPIDInterface",Msg::kError)<<"Will abort here..."<<endl;
00437     assert(false);
00438   }
00439   
00440   //open input file stream
00441   ifstream file(sFileName.c_str());
00442   
00443   //check if file exists
00444   if(!file){
00445     //file does not exist
00446     MSG("NuPIDInterface",Msg::kError)
00447       <<endl<<endl
00448       <<"GetFileNameRoID: file does not exist, filename="
00449       <<sFileName<<endl<<endl;
00450     MSG("NuPIDInterface",Msg::kError)
00451       <<endl<<endl
00452       <<"The RO PID input files are no longer distributed so they"
00453       <<endl
00454       <<"have to be copied to the local site. The NuMuBar/data"
00455       <<endl
00456       <<"directory is the place to copy or soft link the files to"
00457       <<endl
00458       <<"The Cedar/Daikon files are documented here:"
00459       <<endl
00460       <<"    http://home.fnal.gov/~rustem/PhysicsNtuple/"
00461       <<endl
00462       <<"and possibly located here:"
00463       <<endl
00464       <<"    /afs/fnal.gov/files/home/room1/rustem/data/muon-knn"
00465       <<endl<<endl;
00466     MSG("NuPIDInterface",Msg::kError)
00467       <<endl<<endl<<"Asserting here..."<<endl<<endl;
00468     assert(false);
00469   }
00470   else {
00471     //file exists, test is complete so close ifstream
00472     MSG("NuPIDInterface",Msg::kDebug)
00473       <<"File exists with name="<<sFileName<<endl;
00474     file.close();
00475   }
00476   
00477   return sFileName;
00478 }
00479 
00480 //......................................................................
00481 
00482 std::string NuPIDInterface::GetFileNameJeID(const NuEvent& nu) const
00483 {
00484   //get an instance of the code library
00485   const NuLibrary& lib=NuLibrary::Instance();
00486   
00487   string sBase=lib.general.GetReleaseDirToUse("NuMuBar");
00488   string sFileName=sBase+"/NuMuBar/data/pid_evans_cedar_daikon_LE185.root";
00489   Bool_t badVersions=false;
00490   
00491   if (nu.detector==Detector::kNear) {
00492     if (ReleaseType::IsCedar(nu.recoVersion) &&
00493         ReleaseType::IsDaikon(nu.mcVersion)) {
00494       sFileName=sBase+"/NuMuBar/data/pid_evans_cedar_daikon_LE185.root";
00495     }
00496     else {
00497       badVersions=true;
00498     }
00499   }
00500   else if (nu.detector==Detector::kFar) {
00501     badVersions=true;
00502   }
00503   else {
00504     cout<<"Ahhh, detector="<<nu.detector<<" not recognised"<<endl;
00505   }
00506 
00507   if (badVersions) {
00508     MSG("NuPIDInterface",Msg::kError)
00509       <<"Ahhh, ReleaseType not recognised:"<<endl;
00510     this->PrintReleaseTypeEtc(nu);
00511     MSG("NuPIDInterface",Msg::kError)<<"Will abort here..."<<endl;
00512     assert(false);
00513   }
00514 
00515   return sFileName;
00516 }
00517 
00518 //......................................................................
00519 // PoID Officially depracated
00520 //void NuPIDInterface::InitialisePoID(const NuEvent& nu)
00521 //{
00522 //  MSG("NuPIDInterface",Msg::kInfo)
00523 //    <<"InitialisePoID..."<<endl;
00524 //
00525 //  TDirectory* tmpd = gDirectory;
00526 //  MSG("NuPIDInterface",Msg::kDebug)
00527 //    <<"TDirectory before InitialisePoID is:"<<endl;
00528 //  //tmpd->pwd();
00529 //
00530 //  string sFileName=this->GetFileNamePoID(nu);
00531 //  MSG("NuPIDInterface",Msg::kInfo)
00532 //    <<"NuBarPID::LoadPDF with:"<<endl
00533 //    <<" sFileName="<<sFileName<<endl;
00534 //  this->PrintReleaseTypeEtc(nu);
00535 //
00536 //  string sFileNameKin=this->GetFileNamePoIDKin(nu);
00537 //  MSG("NuPIDInterface",Msg::kInfo)
00538 //    <<"NuBarPID::LoadPDF (Kin) with:"<<endl
00539 //    <<" sFileNameKin="<<sFileNameKin<<endl;
00540 //  this->PrintReleaseTypeEtc(nu);
00541 //
00542 //  //load the files
00543 //  poLoadPDF=new LoadPDF(sFileName.c_str());
00544 //  poLoadPDFKin=new LoadPDF(sFileNameKin.c_str());
00545 //
00546 //  //create the calculator object
00547 //  string sCalcType="QYC";
00548 //  MSG("NuPIDInterface",Msg::kInfo)
00549 //    <<"Creating NuBarPID Calculator with sCalcType="<<sCalcType<<endl;
00550 //  poPIDCalculator=new PIDCalculator(sCalcType.c_str()); 
00551 //  
00552 //  //poPIDCalculator=new PIDCalculator("QYC"); 
00553 //
00554 //  //create the second calculator object for kinematics-only PID
00555 //  string sCalcTypeKin="PYC";
00556 //  MSG("NuPIDInterface",Msg::kInfo)
00557 //    <<"Creating NuBarPID kinematics-only calculator with sCalcTypeKin="<<sCalcTypeKin<<endl;
00558 //  poPIDCalculatorKin=new PIDCalculator(sCalcTypeKin.c_str()); 
00559 //
00560 //  MSG("NuPIDInterface",Msg::kDebug)
00561 //    <<"TDirectory after InitialisePoID is:"<<endl;
00562 //  //gDirectory->pwd();
00563 //  MSG("NuPIDInterface",Msg::kDebug)
00564 //    <<"After reseting gDirectory its location is:"<<endl;
00565 //  gDirectory=tmpd;
00566 //  //gDirectory->pwd();
00567 //}
00568 
00569 //......................................................................
00570 
00571 void NuPIDInterface::InitialiseAbID(const NuEvent& nu) 
00572 {
00573   MSG("NuPIDInterface",Msg::kInfo)
00574     <<"InitialiseAbID..."<<endl;
00575 
00576   TDirectory* tmpd = gDirectory;
00577   MSG("NuPIDInterface",Msg::kDebug)
00578     <<"TDirectory before InitialiseAbID is:"<<endl;
00579   //tmpd->pwd();
00580 
00581   string sFileName=this->GetFileNameAbID(nu);
00582   MSG("NuPIDInterface",Msg::kInfo)
00583     <<"InitialiseAbID with:"<<endl
00584     <<" sFileName="<<sFileName<<endl;
00585   this->PrintReleaseTypeEtc(nu);
00586 
00587   //construct the object
00588   abID=new MadAbID();
00589 
00590   //open the file
00591   abID->ReadPDFs(sFileName.c_str());
00592   
00593   MSG("NuPIDInterface",Msg::kDebug)
00594     <<"TDirectory after InitialiseAbID is:"<<endl;
00595   //gDirectory->pwd();
00596   MSG("NuPIDInterface",Msg::kDebug)
00597     <<"After reseting gDirectory its location is:"<<endl;
00598   gDirectory=tmpd;
00599   //gDirectory->pwd();
00600 }
00601 
00602 //......................................................................
00603 
00604 void NuPIDInterface::InitialiseDpID(const NuEvent& nu) 
00605 {
00606   MSG("NuPIDInterface",Msg::kInfo)
00607     <<"InitialiseDpID..."<<endl;
00608 
00609   TDirectory* tmpd = gDirectory;
00610   MSG("NuPIDInterface",Msg::kDebug)
00611     <<"TDirectory before InitialiseDpID is:"<<endl;
00612   //tmpd->pwd();
00613 
00614   //construct the interface to mad
00615   //create these on the heap so that they don't destruct
00616   // - on purpose to avoid a segv on destruction... nice coding!
00617   mad=new NuMadAnalysis();
00618   
00619   //construct the object
00620   dpID=new MadDpID();
00621   
00622   BeamType::BeamType_t btype = static_cast<BeamType::BeamType_t>(nu.beamType);
00623   MAXMSG("NuPIDInterface",Msg::kInfo,1) << "Initializing DpID for btype = " << BeamType::AsString(btype) 
00624   << " (" << btype << ") and " << endl;
00625   this->PrintReleaseTypeEtc(nu);
00626 
00627   // Check for and change RHC to FHC
00628   if (btype == BeamType::kL010z185i_rev) {
00629     btype = BeamType::kL010z185i;
00630   }
00631   MAXMSG("NuPIDInterface",Msg::kInfo,1) << "Initializing DpID for btype = " << BeamType::AsString(btype) 
00632   << " (" << btype << ") and " << endl;
00633     this->PrintReleaseTypeEtc(nu);
00634   
00635   //select the correct pdf
00636   // Force release type to cedar since those are the only training files available. 
00637   int rtype = nu.releaseType;
00638   rtype -= ReleaseType::GetRecoInfo(rtype); // remove the current recoversion
00639   rtype += ReleaseType::kCedar;             // make this appear to be cedar
00640     
00641   dpID->ChoosePDFs(static_cast<Detector::Detector_t>(nu.detector),
00642                    btype,
00643                    ReleaseType::AsString(rtype),
00644                    ReleaseType::AsString(nu.mcVersion));
00645 
00646   MSG("NuPIDInterface",Msg::kDebug)
00647     <<"TDirectory after InitialiseDpID is:"<<endl;
00648   //gDirectory->pwd();
00649   MSG("NuPIDInterface",Msg::kDebug)
00650     <<"After reseting gDirectory its location is:"<<endl;
00651   gDirectory=tmpd;
00652   //gDirectory->pwd();
00653 }
00654 
00655 //......................................................................
00656 
00657 void NuPIDInterface::InitialiseRoID2007(const NuEvent& nu) 
00658 {
00659   MSG("NuPIDInterface",Msg::kInfo)
00660     <<"InitialiseRoID..."<<endl;
00661 
00662   TDirectory* tmpd = gDirectory;
00663   MSG("NuPIDInterface",Msg::kDebug)
00664     <<"TDirectory before InitialiseRoID is:"<<endl;
00665   //tmpd->pwd();
00666 
00667   string sFileName=this->GetFileNameRoID(nu);
00668   MSG("NuPIDInterface",Msg::kInfo)
00669     <<"InitialiseRoID with:"<<endl
00670     <<" sFileName="<<sFileName<<endl;
00671   this->PrintReleaseTypeEtc(nu);
00672   
00673   //get the interface
00674   roInterface=new Anp::Interface();
00675 
00676   //get the registry
00677   roRegistry=new Registry();
00678 
00679   roRegistry->UnLockValues();
00680   roRegistry->Set("FillkNNKNeighbor", int(41));
00681   roRegistry->Set("FillkNNKeyList", "7001, 7010, 7020, 7040");
00682   roRegistry->Set("FillkNNMinViewNPlane", int(5));
00683   //roRegistry->Set("FillkNNFilePath", "train.root");
00684   roRegistry->Set("FillkNNFilePath", sFileName.c_str());
00685   roRegistry->Set("FillkNNKeyBase", int(5200));
00686   roRegistry->Set("QuietInterface", "yes");
00687   roRegistry->Set("FillMuonIdStudy","no");
00688   roRegistry->Set("RunAlgEventErase","no");
00689   
00690   roRegistry->Set("AlgStoreList", "FillHeader FillShower FillTrack FillStrip FillEvent");
00691   roRegistry->Set("AlgSnarlList", "FillMuonId FillTrackGeom FillkNN RunAlgEvent");
00692   roRegistry->Set("AlgEventList","SelectAntiNeutrino");
00693   roRegistry ->Set("PrintConfig","yes");
00694   roRegistry->LockValues();
00695 
00696   // configure RoID interface with registry values
00697   roInterface->Config(*roRegistry);
00698 
00699   MSG("NuPIDInterface",Msg::kDebug)
00700     <<"TDirectory after InitialiseRoID is:"<<endl;
00701   //gDirectory->pwd();
00702   MSG("NuPIDInterface",Msg::kDebug)
00703     <<"After reseting gDirectory its location is:"<<endl;
00704   gDirectory=tmpd;
00705   //gDirectory->pwd();
00706 }
00707 
00708 //......................................................................
00709 
00710 void NuPIDInterface::InitialiseRoID(const NuEvent& nu) 
00711 {
00712   MSG("NuPIDInterface",Msg::kInfo)
00713     <<"InitialiseRoID..."<<endl;
00714 
00715   TDirectory* tmpd = gDirectory;
00716   MSG("NuPIDInterface",Msg::kInfo)
00717     <<"TDirectory before InitialiseRoID is:"<<endl;
00718   tmpd->pwd();
00719 
00720   const string sFileName=this->GetFileNameRoID(nu);
00721   MSG("NuPIDInterface",Msg::kInfo)
00722     <<"InitialiseRoID with:"<<endl
00723     <<" sFileName="<<sFileName<<endl;
00724   this->PrintReleaseTypeEtc(nu);
00725   
00726   //get the interface
00727   roInterface=new Anp::Interface();
00728 
00729   //get an instance of the code library
00730   const NuLibrary& lib=NuLibrary::Instance();
00731   
00732   //get the path to the config file 
00733   const string sBase=lib.general.GetReleaseDirToUse("PhysicsNtuple");                       
00734   const string sConfigFileName=sBase+ "/PhysicsNtuple/Config/Config2008Test.txt";
00735   
00736   //get the registry
00737   roRegistry=new Registry(false);
00738   roRegistry->Set("InterfaceConfigPath",sConfigFileName.c_str());
00739   roRegistry->Set("FillkNNFilePath",sFileName.c_str());
00740   roRegistry ->Set("PrintConfig","yes");
00741   cout<<" inside initialise roInterace"<<endl;
00742   // configure RoID interface with registry values
00743   roInterface->Config(*roRegistry);
00744 
00745   MSG("NuPIDInterface",Msg::kInfo)
00746     <<"TDirectory after InitialiseRoID is:"<<endl;
00747   gDirectory->pwd();
00748   MSG("NuPIDInterface",Msg::kInfo)
00749     <<"After reseting gDirectory its location is:"<<endl;
00750   gDirectory=tmpd;
00751   gDirectory->pwd();
00752 }
00753 //-------------------------------------------------------------------------------------
00754 void NuPIDInterface::InitialiseJmID(const NuEvent& nu)
00755 {
00756   MSG("NuPIDInterface",Msg::kInfo)
00757     <<"InitialiseRoID..."<<endl;
00758 
00759   TDirectory* tmpd = gDirectory;
00760   MSG("NuPIDInterface",Msg::kInfo)
00761     <<"TDirectory before InitialiseJmID is:"<<endl;
00762   tmpd->pwd();
00763 
00764   const string sFileName=this->GetFileNameJmID(nu);
00765   MSG("NuPIDInterface",Msg::kInfo)
00766     <<"InitialiseJmID with:"<<endl
00767     <<" sFileName="<<sFileName<<endl;
00768   this->PrintReleaseTypeEtc(nu);
00769 
00770   //get the interface                                                                                                                                                                                          
00771   jmInterface=new Anp::Interface();
00772 
00773   //get an instance of the code library                                                                                                                                                                        
00774   const NuLibrary& lib=NuLibrary::Instance();
00775 
00776   //get the path to the config file          
00777   const string sBase=lib.general.GetReleaseDirToUse("PhysicsNtuple");                       
00778   const string sConfigFileName=sBase+ "/PhysicsNtuple/Config/ShortConfig2009.txt";
00779   //  const string sConfigFileName="ShortConfig2009.txt";
00780 
00781   //get the registry 
00782   cout<<" inside initialise jmInterface "<<sConfigFileName<<endl;
00783                          
00784   jmRegistry=new Registry(false);
00785   string knn_name2;                                                                                                                      
00786   const string sBaseE=lib.general.GetReleaseDirToUse("NuMuBar");
00787   if(!(sFileName.find("near") ==string::npos))
00788     knn_name2 = "/NuMuBar/data/knn.physics.notrack.near.daikon_04.dogwood1.L010z185i.root";
00789   else if(!(sFileName.find("far") ==string::npos))
00790     knn_name2 = "/NuMuBar/data/knn.physics.notrack.far.daikon_04.dogwood1.L010z185i.root";
00791 
00792   const string knn2 = sBaseE+ knn_name2;
00793 
00794   jmRegistry->Set("EventkNN",*SubReg_FillkNN(99,19400,19000, "19008 19007 19004",knn2,"no"));
00795 
00796   jmRegistry->Set("InterfaceConfigPath",sConfigFileName.c_str());
00797   jmRegistry->Set("FillkNNFilePath",sFileName.c_str());
00798   jmRegistry ->Set("PrintConfig","yes");
00799 
00800   // configure RoID interface with registry values                                                                                                                                                             
00801   jmInterface->Config(*jmRegistry);
00802 
00803   MSG("NuPIDInterface",Msg::kInfo)
00804     <<"TDirectory after InitialiseJmID is:"<<endl;
00805   gDirectory->pwd();
00806   MSG("NuPIDInterface",Msg::kInfo)
00807     <<"After reseting gDirectory its location is:"<<endl;
00808   gDirectory=tmpd;
00809   gDirectory->pwd();
00810 }
00811 
00812 //......................................................................
00813 
00814 void NuPIDInterface::InitialiseRoIDNuMuBar(const NuEvent& nu) 
00815 {
00816   MSG("NuPIDInterface",Msg::kInfo)
00817     <<"InitialiseRoIDNuMuBar..."<<endl;
00818 
00819   TDirectory* tmpd = gDirectory;
00820   MSG("NuPIDInterface",Msg::kDebug)
00821     <<"TDirectory before InitialiseRoIDNuMuBar is:"<<endl;
00822   //tmpd->pwd();
00823 
00824   //string sFileName=this->GetFileNameRoID(nu);
00825   //MSG("NuPIDInterface",Msg::kInfo)
00826   //<<"InitialiseRoIDNuMuBar with:"<<endl
00827   //<<" sFileName="<<sFileName<<endl;
00828   this->PrintReleaseTypeEtc(nu);
00829   
00830   //get the interface
00831   roInterface=new Anp::Interface();
00832 
00833   //get the registry
00834   roRegistry=new Registry();
00835   
00836   roRegistry->UnLockValues();
00837   roRegistry->Set("FillMuonIdStudy","no");
00838   roRegistry->Set("QuietInterface","yes");
00839   roRegistry->Set("RunAlgEventErase","no");
00840   roRegistry->Set("AlgStoreList","FillHeader FillShower FillTrack FillStrip FillEvent");
00841   roRegistry->Set("AlgSnarlList","FillMuonId FillTrackGeom RunAlgEvent");
00842   roRegistry->Set("AlgEventList","SelectAntiNeutrino");
00843   roRegistry ->Set("PrintConfig","yes");
00844   roRegistry->LockValues();
00845   
00846   // configure RoID interface with registry values
00847   roInterface->Config(*roRegistry);
00848   
00849   MSG("NuPIDInterface",Msg::kDebug)
00850     <<"TDirectory after InitialiseRoIDNuMuBar is:"<<endl;
00851   //gDirectory->pwd();
00852   MSG("NuPIDInterface",Msg::kDebug)
00853     <<"After reseting gDirectory its location is:"<<endl;
00854   gDirectory=tmpd;
00855   //gDirectory->pwd();
00856 }
00857 
00858 //......................................................................
00859 
00860 void NuPIDInterface::InitialiseJeID(const NuEvent& nu) 
00861 {
00862   MSG("NuPIDInterface",Msg::kInfo)
00863     <<"InitialiseJeID..."<<endl;
00864 
00865   TDirectory* tmpd = gDirectory;
00866   MSG("NuPIDInterface",Msg::kDebug)
00867     <<"TDirectory before InitialiseJeID is:"<<endl;
00868   //tmpd->pwd();
00869 
00870   string sFileName=this->GetFileNameJeID(nu);
00871   MSG("NuPIDInterface",Msg::kInfo)
00872     <<"InitialiseJeID with:"<<endl
00873     <<" sFileName="<<sFileName<<endl
00874     <<" recoVersion="<<nu.recoVersion
00875     <<", mcVersion="<<nu.mcVersion<<endl;
00876 
00877   //nothing here yet...
00878 
00879   MSG("NuPIDInterface",Msg::kDebug)
00880     <<"TDirectory after InitialiseJeID is:"<<endl;
00881   //gDirectory->pwd();
00882   MSG("NuPIDInterface",Msg::kDebug)
00883     <<"After reseting gDirectory its location is:"<<endl;
00884   gDirectory=tmpd;
00885   //gDirectory->pwd();
00886 }
00887 
00888 //......................................................................
00889 // PoID Officially depracated
00890 //void NuPIDInterface::GetPoID(NuEvent& nu)
00891 //{
00892 //  if (!poLoadPDF) this->InitialisePoID(nu);
00893 //  if (!poLoadPDFKin) this->InitialisePoID(nu);
00894 //
00895 //  //poPIDCalculator->GetNuBarPID(track_q/p*1/track_sigma_q/p,track_nplanes,reco_y,reco_dcosz);
00896 //  
00897 //  //define trk length to be the same as the mad "trklength" variable
00898 //  Int_t trkLength=nu.planeTrkEnd-nu.planeTrkBeg;
00899 //  
00900 //  nu.poID=poPIDCalculator->GetNuBarPID(1./nu.sigqp_qp,trkLength,
00901 //                                       nu.y,nu.trkvtxdcosz,poLoadPDF);
00902 //  //used to use this up to 20070816
00903 //  //nu.poID=poPIDCalculator->GetNuBarPID(1./nu.sigqp_qp,nu.trknplane,
00904 //  //                         nu.y,nu.trkvtxdcosz,poLoadPDF);
00905 //  
00906 //  //(q/p)/(sigma q/p)
00907 //  //Number of planes of track
00908 //  //reconstructed y
00909 //  //reconstructed cosine of angle with respect to beam.
00910 //  
00911 //  //1) that the event is inside the fiducial volume (is_fid==1 in MadMKAnalysis pans)
00912 //  
00913 //  //2) Track quality cut (pitt_trk_qual==1 in MadMKAnalysis pans)
00914 //  
00915 //  //3) Fit probability cut (TMath::Prob(trkchi2,trkndf)>0.1 in MadMKAnalysis pans). 
00916 //  
00917 //  //Kinematics-only pid
00918 //  nu.poIDKin=poPIDCalculatorKin->GetNuBarPID(1./nu.sigqp_qp,trkLength,
00919 //                                             nu.y,nu.trkvtxdcosz,poLoadPDFKin);
00920 //}
00921 
00922 //......................................................................
00923 
00924 void NuPIDInterface::GetAbID(const NtpStRecord& ntp,
00925                            const NtpSREvent& evt,NuEvent& nu)
00926 {
00927   //initialise if required
00928   if (!abID) this->InitialiseAbID(nu);
00929 
00930   //get the value
00931   nu.abID=abID->CalcPID(&evt,&ntp);
00932 }
00933 
00934 //......................................................................
00935 
00936 void NuPIDInterface::GetRoID2007(const NtpStRecord& ntp,
00937                                  const NtpSREvent& evt,NuEvent& nu)
00938 {
00939   //initialise if required
00940   if (!roInterface) this->InitialiseRoID(nu);
00941   
00942   //ro_pass keeps the same value unless it's a new snarl
00943   static Bool_t ro_pass=false;
00944   static Int_t lastSnarl=-999;
00945   MAXMSG("NuPIDInterface",Msg::kDebug,200)
00946     <<"lastSnarl="<<lastSnarl<<", nu.entry="<<nu.entry<<endl;
00947   if (nu.entry!=lastSnarl) {
00948     MAXMSG("NuPIDInterface",Msg::kDebug,200)
00949       <<"NuPIDInterface::GetRoID FillSnarl..."<<endl;
00950     ro_pass=roInterface->FillSnarl(const_cast<NtpStRecord*>(&ntp));
00951   }
00952   lastSnarl=nu.entry;
00953 
00954   Float_t roID=-999;
00955   Float_t roIDNuMuBar=0;//if not NuMuBar then no key in map below
00956   Float_t relativeAngle=-999;
00957 
00958   //variables to record when a variable has been found
00959   Bool_t got5240=false;
00960   Bool_t got5524=false;
00961   Bool_t got5533=false;
00962 
00963   if (ro_pass) {    
00964     map<short, float> romap=roInterface->GetData
00965       (const_cast<NtpSREvent*>(&evt));    
00966     
00967     //Float_t roid_trkplanefrac; // fraction of track planes/total planes
00968     //Float_t roid_lowphhiph;    // fraction of "lower ph"/"higher ph"
00969     //Float_t roid_trkphfracend; // fraction of track ph at end/track ph
00970     //Float_t roid_trkphfracclose; // fraction of track ph/ph of strips
00971     // around track
00972 
00973     for (map<short,float>::const_iterator rit=romap.begin();
00974          rit!=romap.end(); ++rit) {      
00975       Int_t key=0;
00976       Float_t val=0;
00977       key=rit->first;
00978       val=rit->second;
00979       
00980       MAXMSG("NuPIDInterface",Msg::kDebug,1000)
00981         <<"key="<<key<<", val="<<val<<endl;
00982 
00983       if (key==5240) {
00984         roID=val;
00985         got5240=true;
00986       }
00987 
00988       if (key==5524) {
00989         roIDNuMuBar=val;
00990         got5524=true;
00991       }
00992       
00993       if (key==5533) {
00994         relativeAngle=val;
00995         got5533=true;
00996       }
00997       
00998       if (got5240 && got5524 && got5533) {
00999         MAXMSG("NuPIDInterface",Msg::kDebug,50)
01000           <<"Jumping out of loop over ROID map..."<<endl;
01001         break;//jump out of loop once all are found
01002       }
01003       
01004       //if (key==7000) roid_trkplanefrac=val;
01005       //if (key==7020) roid_lowphhiph=val;
01006       //if (key==7030) roid_trkphfracend=val;
01007       //if (key==7040) roid_trkphfracclose=val;
01008     }
01009   }
01010   else {
01011     MAXMSG("NuPIDInterface",Msg::kInfo,100)
01012       <<"Event did not pass RO ID FillSnarl"<<endl;
01013   }
01014 
01015   //store the variables
01016   nu.roID=roID;
01017   nu.roIDNuMuBar=roIDNuMuBar;
01018   nu.relativeAngle=relativeAngle;
01019 
01020   MAXMSG("NuPIDInterface",Msg::kDebug,100)
01021     <<"GetRoID: roID="<<nu.roID
01022     <<", roIDNuMuBar"<<nu.roIDNuMuBar
01023     <<", relativeAngle"<<nu.relativeAngle
01024     <<endl;
01025 }
01026 
01027 //......................................................................
01028 
01029 void NuPIDInterface::GetRoID(const NtpStRecord& ntp,
01030                              const NtpSREvent& evt,NuEvent& nu)
01031 {
01032   //initialise if required
01033   if (!roInterface) this->InitialiseRoID(nu);
01034 
01035   //variable to turn on all the useful messages if required
01036   Msg::LogLevel_t logLevel=Msg::kDebug;
01037 
01038   //have to protect against PhysicsNtuple altering the global state
01039   TDirectory* tmpd = gDirectory;
01040   MSG("NuPIDInterface",Msg::kDebug)
01041     <<"TDirectory before GetRoID is:"<<endl;
01042   //tmpd->pwd();
01043 
01044   //ro_pass keeps the same value unless it's a new snarl
01045   static Bool_t ro_pass=false;
01046   static Int_t lastSnarl=-999;
01047   MAXMSG("NuPIDInterface",logLevel,100)
01048     <<"lastSnarl="<<lastSnarl<<", nu.entry="<<nu.entry<<endl;
01049   if (nu.entry!=lastSnarl) {
01050     MAXMSG("NuPIDInterface",logLevel,50)
01051       <<"NuPIDInterface::GetRoID FillSnarl..."<<endl;
01052     ro_pass=roInterface->FillSnarl(const_cast<NtpStRecord*>(&ntp));
01053   }
01054   lastSnarl=nu.entry;
01055 
01056   //get the variables
01057   if (ro_pass) {    
01058     //get the roid for the evt, trk1, trk2 and trk3
01059     this->GetRoIDEvt(evt,nu);
01060     this->GetRoIDPrimaryTrk(ntp,evt,nu);
01061     this->GetRoIDSecondTrk(ntp,evt,nu);
01062     this->GetRoIDThirdTrk(ntp,evt,nu);
01063   }
01064   else {
01065     MAXMSG("NuPIDInterface",Msg::kWarning,100)
01066       <<"Event did not pass RO ID FillSnarl"<<endl;
01067   }
01068 
01069   //if the event is not identified as a NuMuBar then 
01070   //set the value of roIDNuMuBar to be 0 (rather than just the default)
01071   //check that it isn't -999 because then there was no track
01072   if (nu.roIDNuMuBarEvt<0 && 
01073       nu.roIDNuMuBarEvt!=-999) nu.roIDNuMuBarEvt=0;
01074   if (nu.roIDNuMuBar1<0 && nu.roIDNuMuBar1!=-999) nu.roIDNuMuBar1=0;
01075   if (nu.roIDNuMuBar2<0 && nu.roIDNuMuBar2!=-999) nu.roIDNuMuBar2=0;
01076   if (nu.roIDNuMuBar3<0 && nu.roIDNuMuBar3!=-999) nu.roIDNuMuBar3=0;
01077   
01078   //copy the variables across to the "best" roid variables
01079   //these variables allow a level of indirection and can be 
01080   //changed to the "best" value at a later stage
01081   //for now, just use the evt
01082 
01083   nu.roID=nu.roIDEvt;
01084   nu.knn01TrkActivePlanes=nu.knn01TrkActivePlanesEvt;
01085   nu.knn10TrkMeanPh=nu.knn10TrkMeanPhEvt;
01086   nu.knn20LowHighPh=nu.knn20LowHighPhEvt;
01087   nu.knn40TrkPhFrac=nu.knn40TrkPhFracEvt;
01088   nu.roIDNuMuBar=nu.roIDNuMuBarEvt;
01089   nu.relativeAngle=nu.relativeAngleEvt;
01090 
01091   MAXMSG("NuPIDInterface",logLevel,100)
01092     <<"GetRoIDEvt: roID="<<nu.roIDEvt
01093     <<", roIDNuMuBar="<<nu.roIDNuMuBarEvt
01094     <<", relativeAngle="<<nu.relativeAngleEvt
01095     <<endl
01096     <<"   knn: TrkActivePlanes="<<nu.knn01TrkActivePlanesEvt
01097     <<", TrkMeanPh="<<nu.knn10TrkMeanPhEvt
01098     <<", LowHighPh="<<nu.knn20LowHighPhEvt
01099     <<", TrkPhFrac="<<nu.knn40TrkPhFracEvt
01100     <<endl;
01101   MAXMSG("NuPIDInterface",logLevel,100)
01102     <<"GetRoIDTrk1: roID="<<nu.roID1
01103     <<", roIDNuMuBar="<<nu.roIDNuMuBar1
01104     <<", relativeAngle="<<nu.relativeAngle1
01105     <<endl
01106     <<"   knn: TrkActivePlanes="<<nu.knn01TrkActivePlanes1
01107     <<", TrkMeanPh="<<nu.knn10TrkMeanPh1
01108     <<", LowHighPh="<<nu.knn20LowHighPh1
01109     <<", TrkPhFrac="<<nu.knn40TrkPhFrac1
01110     <<endl;
01111   MAXMSG("NuPIDInterface",logLevel,100)
01112     <<"GetRoIDTrk2: roID="<<nu.roID2
01113     <<", roIDNuMuBar="<<nu.roIDNuMuBar2
01114     <<", relativeAngle="<<nu.relativeAngle2
01115     <<endl
01116     <<"   knn: TrkActivePlanes="<<nu.knn01TrkActivePlanes2
01117     <<", TrkMeanPh="<<nu.knn10TrkMeanPh2
01118     <<", LowHighPh="<<nu.knn20LowHighPh2
01119     <<", TrkPhFrac="<<nu.knn40TrkPhFrac2
01120     <<endl;
01121   MAXMSG("NuPIDInterface",logLevel,100)
01122     <<"GetRoIDTrk3: roID="<<nu.roID3
01123     <<", roIDNuMuBar="<<nu.roIDNuMuBar3
01124     <<", relativeAngle="<<nu.relativeAngle3
01125     <<endl
01126     <<"   knn: TrkActivePlanes="<<nu.knn01TrkActivePlanes3
01127     <<", TrkMeanPh="<<nu.knn10TrkMeanPh3
01128     <<", LowHighPh="<<nu.knn20LowHighPh3
01129     <<", TrkPhFrac="<<nu.knn40TrkPhFrac3
01130     <<endl<<endl;
01131 
01132   MSG("NuPIDInterface",Msg::kDebug)
01133     <<endl<<endl<<endl<<"TDirectory after GetRoID is:"<<endl;
01134   //gDirectory->pwd();
01135   MSG("NuPIDInterface",Msg::kDebug)
01136     <<endl<<endl<<endl<<"After reseting gDirectory it is:"<<endl;
01137   gDirectory=tmpd;
01138   //gDirectory->pwd();
01139 }
01140 
01141 //......................................................................
01142 
01143 void NuPIDInterface::GetRoIDEvt(const NtpSREvent& evt,
01144                                 NuEvent& nu) const
01145 {
01146   if (!roInterface) {
01147     MSG("NuPIDInterface",Msg::kError)
01148       <<"NuPIDInterface: roInterface not initialised = "<<roInterface
01149       <<endl<<"Will assert here..."<<endl;
01150     assert(false);
01151   }
01152   //NtpStRecord* st=const_cast<NtpStRecord*>(&ntp);//nasty!
01153   NtpSREvent* pevt=const_cast<NtpSREvent*>(&evt);//nasty!
01154   
01155   //float GetVar(const std::string &, TObject *) const; 
01156   nu.roIDEvt=roInterface->GetVar("knn_pid",pevt);
01157   nu.knn01TrkActivePlanesEvt=roInterface->GetVar("knn_01",pevt);
01158   nu.knn10TrkMeanPhEvt=roInterface->GetVar("knn_10",pevt);
01159   nu.knn20LowHighPhEvt=roInterface->GetVar("knn_20",pevt);
01160   nu.knn40TrkPhFracEvt=roInterface->GetVar("knn_40",pevt);
01161   nu.roIDNuMuBarEvt=roInterface->GetVar("numubar",pevt);
01162   nu.relativeAngleEvt=roInterface->GetVar("rel_ang",pevt);
01163 
01164 }
01165 
01166 //......................................................................
01167 
01168 void NuPIDInterface::GetRoIDPrimaryTrk(const NtpStRecord& ntp,
01169                                        const NtpSREvent& evt,
01170                                        NuEvent& nu) const
01171 {
01172   if (!roInterface) {
01173     MSG("NuPIDInterface",Msg::kError)
01174       <<"NuPIDInterface: roInterface not initialised = "<<roInterface
01175       <<endl<<"Will assert here..."<<endl;
01176     assert(false);
01177   }
01178   //get an instance of the code library
01179   const NuLibrary& lib=NuLibrary::Instance();
01180   
01181   //get pointer to primary trk
01182   const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,0);
01183   
01184   //check if track exists
01185   if (ptrk) {
01186     //have to const cast to pass to GetVar
01187     NtpSRTrack* ptrkb=const_cast<NtpSRTrack*>(ptrk);//nasty!
01188     
01189     nu.roID1=roInterface->GetVar("knn_pid",ptrkb);
01190     nu.knn01TrkActivePlanes1=roInterface->GetVar("knn_01",ptrkb);
01191     nu.knn10TrkMeanPh1=roInterface->GetVar("knn_10",ptrkb);
01192     nu.knn20LowHighPh1=roInterface->GetVar("knn_20",ptrkb);
01193     nu.knn40TrkPhFrac1=roInterface->GetVar("knn_40",ptrkb);
01194     nu.roIDNuMuBar1=roInterface->GetVar("numubar",ptrkb);
01195     nu.relativeAngle1=roInterface->GetVar("rel_ang",ptrkb);
01196   }
01197 }
01198 
01199 //......................................................................
01200 
01201 void NuPIDInterface::GetRoIDSecondTrk(const NtpStRecord& ntp,
01202                                       const NtpSREvent& evt,
01203                                       NuEvent& nu) const
01204 {
01205   if (!roInterface) {
01206     MSG("NuPIDInterface",Msg::kError)
01207       <<"NuPIDInterface: roInterface not initialised = "<<roInterface
01208       <<endl<<"Will assert here..."<<endl;
01209     assert(false);
01210   }
01211 
01212   //get an instance of the code library
01213   const NuLibrary& lib=NuLibrary::Instance();
01214   
01215   //get pointer to second trk
01216   const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,1);
01217   
01218   //check if track exists
01219   if (ptrk) {
01220     //have to const cast to pass to GetVar
01221     NtpSRTrack* ptrkb=const_cast<NtpSRTrack*>(ptrk);//nasty!
01222     
01223     nu.roID2=roInterface->GetVar("knn_pid",ptrkb);
01224     nu.knn01TrkActivePlanes2=roInterface->GetVar("knn_01",ptrkb);
01225     nu.knn10TrkMeanPh2=roInterface->GetVar("knn_10",ptrkb);
01226     nu.knn20LowHighPh2=roInterface->GetVar("knn_20",ptrkb);
01227     nu.knn40TrkPhFrac2=roInterface->GetVar("knn_40",ptrkb);
01228     nu.roIDNuMuBar2=roInterface->GetVar("numubar",ptrkb);
01229     nu.relativeAngle2=roInterface->GetVar("rel_ang",ptrkb);
01230   }
01231 }
01232 
01233 //......................................................................
01234 
01235 void NuPIDInterface::GetRoIDThirdTrk(const NtpStRecord& ntp,
01236                                      const NtpSREvent& evt,
01237                                      NuEvent& nu) const
01238 {
01239   if (!roInterface) {
01240     MSG("NuPIDInterface",Msg::kError)
01241       <<"NuPIDInterface: roInterface not initialised = "<<roInterface
01242       <<endl<<"Will assert here..."<<endl;
01243     assert(false);
01244   }
01245 
01246   //get an instance of the code library
01247   const NuLibrary& lib=NuLibrary::Instance();
01248   
01249   //get pointer to third trk
01250   const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,2);
01251   
01252   //check if track exists
01253   if (ptrk) {
01254     //have to const cast to pass to GetVar
01255     NtpSRTrack* ptrkb=const_cast<NtpSRTrack*>(ptrk);//nasty!
01256     
01257     nu.roID3=roInterface->GetVar("knn_pid",ptrkb);
01258     nu.knn01TrkActivePlanes3=roInterface->GetVar("knn_01",ptrkb);
01259     nu.knn10TrkMeanPh3=roInterface->GetVar("knn_10",ptrkb);
01260     nu.knn20LowHighPh3=roInterface->GetVar("knn_20",ptrkb);
01261     nu.knn40TrkPhFrac3=roInterface->GetVar("knn_40",ptrkb);
01262     nu.roIDNuMuBar3=roInterface->GetVar("numubar",ptrkb);
01263     nu.relativeAngle3=roInterface->GetVar("rel_ang",ptrkb);
01264   }
01265 }
01266 
01267 
01268 //......................................................................
01269 
01270 void NuPIDInterface::GetJmID(const NtpStRecord& ntp,
01271                              const NtpSREvent& evt,NuEvent& nu)
01272 {
01273 
01274   //initialise if required
01275   if (!jmInterface) this->InitialiseJmID(nu);
01276 
01277   //variable to turn on all the useful messages if required
01278   Msg::LogLevel_t logLevel=Msg::kDebug;
01279 
01280   //have to protect against PhysicsNtuple altering the global state
01281   TDirectory* tmpd = gDirectory;
01282   MSG("NuPIDInterface",Msg::kDebug)
01283     <<"TDirectory before GetJmID is:"<<endl;
01284   //tmpd->pwd();
01285 
01286   //jm_pass keeps the same value unless it's a new snarl
01287   static Bool_t jm_pass=false;
01288   static Int_t lastSnarl=-999;
01289   MAXMSG("NuPIDInterface",logLevel,100)
01290     <<"lastSnarl="<<lastSnarl<<", nu.entry="<<nu.entry<<endl;
01291   if (nu.entry!=lastSnarl) {
01292     MAXMSG("NuPIDInterface",logLevel,50)
01293       <<"NuPIDInterface::GetJmID FillSnarl..."<<endl;
01294     jm_pass=jmInterface->FillSnarl(const_cast<NtpStRecord*>(&ntp));
01295   }
01296   lastSnarl=nu.entry;
01297   //get the variables
01298   if (jm_pass) {    
01299     //get the roid for the evt, trk1, trk2 and trk3
01300     this->GetJmIDEvt(evt,nu);
01301     this->GetJmIDPrimaryTrk(ntp,evt,nu);
01302     this->GetJmIDSecondTrk(ntp,evt,nu);
01303     this->GetJmIDThirdTrk(ntp,evt,nu);
01304   }
01305   else {
01306     MAXMSG("NuPIDInterface",Msg::kWarning,100)
01307       <<"Event did not pass JM ID FillSnarl"<<endl;
01308   }
01309   
01310   //copy the variables across to the "best" roid variables
01311   //these variables allow a level of indirection and can be 
01312   //changed to the "best" value at a later stage
01313   //for now, just use the evt
01314 
01315   nu.jmID=nu.jmIDEvt;
01316   nu.jmTrackPlane=nu.jmTrackPlaneEvt;
01317   nu.jmMeanPh=nu.jmMeanPhEvt;
01318   nu.jmEndPh=nu.jmEndPhEvt;
01319   nu.jmScatteringU=nu.jmScatteringUEvt;
01320   nu.jmScatteringV=nu.jmScatteringVEvt;
01321   nu.jmScatteringUV=nu.jmScatteringUVEvt;
01322   nu.jmEventknnID=nu.jmEventknnIDEvt;
01323   nu.jmEventknn208=nu.jmEventknn208Evt;
01324   nu.jmEventknn207=nu.jmEventknn207Evt;
01325   nu.jmEventknn206=nu.jmEventknn206Evt;
01326   nu.jmEventknn205=nu.jmEventknn205Evt;
01327   nu.jmEventknn204=nu.jmEventknn204Evt;
01328 
01329   gDirectory=tmpd;
01330   //  MAXMSG("NuPIDInterface",logLevel,100)
01331   //gDirectory->pwd();
01332 }
01333 
01334 //......................................................................
01335 
01336 void NuPIDInterface::GetJmIDEvt(const NtpSREvent& evt,
01337                                 NuEvent& nu) const
01338 {
01339   if (!jmInterface) {
01340     MSG("NuPIDInterface",Msg::kError)
01341       <<"NuPIDInterface: jmInterface not initialised = "<<jmInterface
01342       <<endl<<"Will assert here..."<<endl;
01343     assert(false);
01344   }
01345   //NtpStRecord* st=const_cast<NtpStRecord*>(&ntp);//nasty!
01346   NtpSREvent* pevt=const_cast<NtpSREvent*>(&evt);//nasty!
01347   
01348   //float GetVar(const std::string &, TObject *) const; 
01349   nu.jmIDEvt=     jmInterface->GetVar("shortid",pevt);
01350   nu.jmTrackPlaneEvt= jmInterface->GetVar("ntrkpl" ,pevt);
01351   nu.jmMeanPhEvt= jmInterface->GetVar("meanph" ,pevt);
01352   nu.jmEndPhEvt=  jmInterface->GetVar("endph"  ,pevt);
01353   nu.jmScatteringUEvt=  jmInterface->GetVar("scatu"  ,pevt);
01354   nu.jmScatteringVEvt=  jmInterface->GetVar("scatv"  ,pevt);
01355   nu.jmScatteringUVEvt= jmInterface->GetVar("scatuv" ,pevt);
01356   nu.jmEventknnIDEvt= jmInterface->GetVar("evtid"  ,pevt);
01357   nu.jmEventknn208Evt=jmInterface->GetVar("evt208" ,pevt);
01358   nu.jmEventknn207Evt=jmInterface->GetVar("evt207" ,pevt);
01359   nu.jmEventknn206Evt=jmInterface->GetVar("evt206" ,pevt);
01360   nu.jmEventknn205Evt=jmInterface->GetVar("evt205" ,pevt);
01361   nu.jmEventknn204Evt=jmInterface->GetVar("evt204" ,pevt);
01362 
01363 }
01364 
01365 //......................................................................
01366 
01367 void NuPIDInterface::GetJmIDPrimaryTrk(const NtpStRecord& ntp,
01368                                        const NtpSREvent& evt,
01369                                        NuEvent& nu) const
01370 {
01371   if (!jmInterface) {
01372     MSG("NuPIDInterface",Msg::kError)
01373       <<"NuPIDInterface: jmInterface not initialised = "<<jmInterface
01374       <<endl<<"Will assert here..."<<endl;
01375     assert(false);
01376   }
01377   //get an instance of the code library
01378   const NuLibrary& lib=NuLibrary::Instance();
01379   
01380   //get pointer to primary trk
01381   const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,0);
01382   
01383   //check if track exists
01384   if (ptrk) {
01385     //have to const cast to pass to GetVar
01386     NtpSRTrack* ptrkb=const_cast<NtpSRTrack*>(ptrk);//nasty!
01387 
01388     nu.jmID1=jmInterface->GetVar("shortid",ptrkb);
01389     nu.jmTrackPlane1=jmInterface->GetVar("ntrkpl",ptrkb);
01390     nu.jmMeanPh1=jmInterface->GetVar("meanph",ptrkb);
01391     nu.jmEndPh1=jmInterface->GetVar("endph",ptrkb);
01392     nu.jmScatteringU1=jmInterface->GetVar("scatu",ptrkb);
01393     nu.jmScatteringV1=jmInterface->GetVar("scatv",ptrkb);
01394   }
01395 }
01396 
01397 //......................................................................
01398 
01399 void NuPIDInterface::GetJmIDSecondTrk(const NtpStRecord& ntp,
01400                                       const NtpSREvent& evt,
01401                                       NuEvent& nu) const
01402 {
01403   if (!jmInterface) {
01404     MSG("NuPIDInterface",Msg::kError)
01405       <<"NuPIDInterface: jmInterface not initialised = "<<jmInterface
01406       <<endl<<"Will assert here..."<<endl;
01407     assert(false);
01408   }
01409 
01410   //get an instance of the code library
01411   const NuLibrary& lib=NuLibrary::Instance();
01412   
01413   //get pointer to second trk
01414   const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,1);
01415   
01416   //check if track exists
01417   if (ptrk) {
01418     //have to const cast to pass to GetVar
01419     NtpSRTrack* ptrkb=const_cast<NtpSRTrack*>(ptrk);//nasty!
01420     
01421     nu.jmID2=jmInterface->GetVar("shortid",ptrkb);
01422     nu.jmTrackPlane2=jmInterface->GetVar("ntrkpl",ptrkb);
01423     nu.jmMeanPh2=jmInterface->GetVar("meanph",ptrkb);
01424     nu.jmEndPh2=jmInterface->GetVar("endph",ptrkb);
01425     nu.jmScatteringU2=jmInterface->GetVar("scatu",ptrkb);
01426     nu.jmScatteringV2=jmInterface->GetVar("scatv",ptrkb);
01427   }
01428 }
01429 
01430 //......................................................................
01431 
01432 void NuPIDInterface::GetJmIDThirdTrk(const NtpStRecord& ntp,
01433                                      const NtpSREvent& evt,
01434                                      NuEvent& nu) const
01435 {
01436   if (!jmInterface) {
01437     MSG("NuPIDInterface",Msg::kError)
01438       <<"NuPIDInterface: jmInterface not initialised = "<<jmInterface
01439       <<endl<<"Will assert here..."<<endl;
01440     assert(false);
01441   }
01442 
01443   //get an instance of the code library
01444   const NuLibrary& lib=NuLibrary::Instance();
01445   
01446   //get pointer to third trk
01447   const NtpSRTrack* ptrk=lib.reco.GetTrackWithIndexX(ntp,evt,2);
01448   
01449   //check if track exists
01450   if (ptrk) {
01451     //have to const cast to pass to GetVar
01452     NtpSRTrack* ptrkb=const_cast<NtpSRTrack*>(ptrk);//nasty!
01453     nu.jmID3=jmInterface->GetVar("shortid",ptrkb);
01454     nu.jmTrackPlane3=jmInterface->GetVar("ntrkpl",ptrkb);
01455     nu.jmMeanPh3=jmInterface->GetVar("meanph",ptrkb);
01456     nu.jmEndPh3=jmInterface->GetVar("endph",ptrkb);
01457     nu.jmScatteringU3=jmInterface->GetVar("scatu",ptrkb);
01458     nu.jmScatteringV3=jmInterface->GetVar("scatv",ptrkb);
01459    
01460   }
01461 }
01462 //--------------------------------------------------------------------------------------------                                                                                                                 
01463 Registry*  NuPIDInterface::SubReg_FillkNN(int k, int keybase, int keytrue, string keylist, string shortknn, string usetrack)
01464 {
01465   Registry *kreg = new Registry(false);
01466 
01467   kreg -> Set("PrintConfig", "yes");
01468   kreg -> Set("AlgSnarlName", "FillkNN");
01469   kreg -> Set("FillkNNKeySignal", int(1));      // for historical reasons muon have type=1                                                                                                               
01470   kreg -> Set("FillkNNKNeighbor", k);
01471   kreg -> Set("FillkNNKNeighborMod", int(5));
01472   kreg -> Set("FillkNNAddEvent", "no");
01473   kreg -> Set("FillkNNUseTrack", usetrack.c_str());
01474   kreg -> Set("FillkNNTrim", "yes");
01475   kreg -> Set("FillkNNTrimDelta", int(1000));
01476   kreg -> Set("FillkNNFilePath", shortknn.c_str());
01477   kreg -> Set("FillkNNKeyTruth", keytrue);
01478   kreg -> Set("FillkNNKeyBase", keybase);
01479   kreg -> Set("FillkNNKeyList", keylist.c_str());
01480   kreg -> Set("FillShortEventKeyPass",14001);
01481   kreg -> Set("FillShortEventKeyPass2",14001); // no access to 11001
01482   kreg -> Set("FillkNN", "no");
01483 
01484   return kreg;
01485 }
01486 
01487 //......................................................................
01488 
01489 void NuPIDInterface::GetRoIDNuMuBar(const NtpStRecord& ntp,
01490                                     const NtpSREvent& evt,NuEvent& nu)
01491 {
01492   //initialise if required
01493   if (!roInterface) this->InitialiseRoIDNuMuBar(nu);
01494   
01495   //ro_pass keeps the same value unless it's a new snarl
01496   static Bool_t ro_pass=false;
01497   static Int_t lastSnarl=-999;
01498   MAXMSG("NuPIDInterface",Msg::kDebug,200)
01499     <<"lastSnarl="<<lastSnarl<<", nu.entry="<<nu.entry<<endl;
01500   if (nu.entry!=lastSnarl) {
01501     MAXMSG("NuPIDInterface",Msg::kDebug,200)
01502       <<"NuPIDInterface::GetRoID FillSnarl..."<<endl;
01503     ro_pass=roInterface->FillSnarl(const_cast<NtpStRecord*>(&ntp));
01504   }
01505   lastSnarl=nu.entry;
01506 
01507   Float_t roIDNuMuBar=0;//if not NuMuBar then no key in map below
01508   Float_t relativeAngle=-999;
01509 
01510   Bool_t got5524=false;
01511   Bool_t got5533=false;
01512 
01513   if (ro_pass) {    
01514     map<short, float> romap=roInterface->GetData
01515       (const_cast<NtpSREvent*>(&evt));    
01516     
01517     MAXMSG("NuPIDInterface",Msg::kDebug,200)
01518       <<"NuPIDInterface: Looping over key/value pairs..."<<endl;
01519 
01520     //loop over the map of data returned
01521     for (map<short,float>::const_iterator rit=romap.begin();
01522          rit!=romap.end(); ++rit) {      
01523       Int_t key=0;
01524       Float_t val=0;
01525       key=rit->first;
01526       val=rit->second;
01527       
01528       MAXMSG("NuPIDInterface",Msg::kDebug,100)
01529         <<"key="<<key<<", val="<<val<<endl;
01530 
01531       if (key==5524) {
01532         roIDNuMuBar=val;
01533         got5524=true;
01534       }
01535       
01536       if (key==5533) {
01537         relativeAngle=val;
01538         got5533=true;
01539       }
01540 
01541       if (got5524 && got5533) {
01542         break;//jump out of loop once all are found
01543       }
01544 
01545     }
01546   }
01547   else {
01548     MAXMSG("NuPIDInterface",Msg::kWarning,100)
01549       <<"Event did not pass RO ID"<<endl;
01550   }
01551   
01552   nu.roIDNuMuBar=roIDNuMuBar;
01553   nu.relativeAngle=relativeAngle;
01554 }
01555 
01556 //......................................................................
01557 
01558 void NuPIDInterface::GetDpID(const NtpStRecord& ntp,
01559                              const NtpSREvent& /*evt*/,NuEvent& nu)
01560 {
01561   //initialise if required
01562   if (!dpID) this->InitialiseDpID(nu);
01563   
01564   //give the ntp to the mad interface
01565   mad->SetEntry(&ntp);
01566   
01567   //get the value
01568   nu.dpID=dpID->CalcPID(mad,nu.evt,0);
01569 }
01570 
01571 //......................................................................
01572 
01573 void NuPIDInterface::PrintReleaseTypeEtc(const NuEvent& nu) const
01574 {
01575   MSG("NuPIDInterface",Msg::kInfo)
01576     <<"ReleaseType="
01577     <<NuUtilities::PrintRelease(nu.releaseType)
01578     <<", detector="<<Detector::AsString
01579     (static_cast<Detector::Detector_t>(nu.detector))<<endl;
01580 }
01581 
01582 //......................................................................
01583