12 #include "Utilities/AssociationUtil.h" 32 #include "TMVA/Factory.h" 33 #include "TMVA/Tools.h" 34 #include "TMVA/Reader.h" 39 #include <boost/algorithm/string.hpp> 40 #include <boost/foreach.hpp> 41 #include <boost/tokenizer.hpp> 44 #define MAX_BUFFER 100000 49 #define FUZZYKVERTEXOUTPUT rb::Prong 126 produces< std::vector<rvp::RVP> >();
127 produces< art::Assns<rvp::RVP,rb::Cluster> >();
168 typedef std::unique_ptr<std::vector<rvp::RVP> > RVPs_t;
169 typedef std::unique_ptr<art::Assns<rvp::RVP,rb::Cluster> > RVPAssn_t;
170 RVPs_t pidcol (
new std::vector<rvp::RVP>);
184 for(
unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
191 const unsigned int nCells = slice.
NCell();
195 double longestTrack = 0.0;
196 int longestTrack_ncells = 0;
198 if (fmTrack.isValid()){
199 std::vector<art::Ptr<rb::Track>> tracks = fmTrack.at(sliceIdx);
201 for(
unsigned int n = 0;
n < tracks.size(); ++
n){
204 if(L > longestTrack){
206 longestTrack_ncells = current_track->
NCell();
212 double recoEnergy = 0.0;
216 const unsigned int ncell_hits = slice.
NCell();
219 std::vector<rb::RecoHit> all_reco_hits;
221 for(
unsigned int icell=0; icell<ncell_hits; ++icell){
227 recoEnergy += rhit.
GeV();
235 if (fmVertex.isValid()){
236 std::vector<art::Ptr<rb::Vertex>> vert = fmVertex.at(sliceIdx);
241 if (vert.size() != 1)
continue;
246 std::vector<double> PlaneSumEnergy;
247 double totalE_road2sigma = 0.0;
250 for(
int iplane=0; iplane<lengthPlanes; ++iplane ){
257 std::vector<double> Xxview;
258 std::vector<double> Exview;
259 std::vector<double> Yyview;
260 std::vector<double> Eyview;
266 for(
unsigned int icell=0; icell<ncell_hits; ++icell){
272 if( Dplane<0 )
std::cout<<
"plane # of this cell is beyond the minimal plane."<<
std::endl;
274 if( Dplane==iplane ){
276 EXcell += rhit.
GeV()*rhit.
X();
277 EYcell += rhit.
GeV()*rhit.
Y();
280 Xxview.push_back(rhit.
X());
281 Exview.push_back(rhit.
GeV());
284 Yyview.push_back(rhit.
Y());
285 Eyview.push_back(rhit.
GeV());
292 cellX = EXcell/Ecell;
293 cellY = EYcell/Ecell;
296 PlaneSumEnergy.push_back(Ecell);
298 const int Nxviews=Xxview.size();
299 const int Nyviews=Yyview.size();
301 for(
int ixv=0; ixv<Nxviews; ++ixv ){
302 if(
fabs(Xxview[ixv]-cellX) < 2.*2. ) totalE_road2sigma+=Exview[ixv];
303 if(
fabs(Xxview[ixv]-cellX) < 2.*3. ) Esig_3+=Exview[ixv];
307 for(
int iyv=0; iyv<Nyviews; ++iyv ){
308 if(
fabs(Yyview[iyv]-cellY) < 2.*2. ) totalE_road2sigma+=Eyview[iyv];
309 if(
fabs(Yyview[iyv]-cellY) < 2.*3. ) Esig_3+=Eyview[iyv];
315 double efrac_2slides = -1.0;
316 double efrac_6slides = -1.0;
317 double Eplane20 = 0.0;
318 int nplanes=PlaneSumEnergy.size();
319 if (recoEnergy > 0.0){
321 if(
i<20 ) Eplane20 += PlaneSumEnergy[
i];
325 const double e_frac2=(PlaneSumEnergy[
i]+PlaneSumEnergy[
i+1])/recoEnergy;
326 if( efrac_2slides<e_frac2 ) efrac_2slides=e_frac2; }
327 else if( nplanes<=2 ) efrac_2slides=1.;
330 const double e_frac6=(PlaneSumEnergy[
i]+PlaneSumEnergy[
i+1]+PlaneSumEnergy[
i+2]+PlaneSumEnergy[
i+3]+PlaneSumEnergy[
i+4]+PlaneSumEnergy[
i+5])/recoEnergy;
331 if( efrac_6slides<e_frac6 ) efrac_6slides=e_frac6; }
332 else if( nplanes<=6 ) efrac_6slides=1.;
343 std::vector<art::Ptr<FUZZYKVERTEXOUTPUT>> prongs = fmProng3D.at(0);
344 std::vector<art::Ptr<FUZZYKVERTEXOUTPUT>> SelectedProng;
345 std::vector<art::Ptr<FUZZYKVERTEXOUTPUT>> SelectedProng2D;
348 SelectedProng = fmProng3D.at(0);
349 SelectedProng2D = fmProng2D.at(0);
352 for(
unsigned int ip=0;
ip<prongs.size(); ++
ip ){
354 if(current_prong->IsNoise() )
continue;
356 SelectedProng.push_back(current_prong);
358 SelectedProng2D.push_back(current_prong);
362 int nprongs = SelectedProng.size();
363 double Ebalance = 1.0;
364 std::vector<double> ProngEnergy;
365 std::vector<double> ProngEnergyXView;
366 std::vector<double> ProngEnergyYView;
369 double EprongX = 0.0;
370 double EprongY = 0.0;
373 const unsigned int nselectedprong_cells = current_prong->NCell();
375 for(
unsigned int icell=0; icell<nselectedprong_cells; ++icell){
377 const rb::RecoHit rhit = current_prong->RecoHit(icell);
380 Eprong += rhit.
GeV();
384 ProngEnergy.push_back(Eprong);
385 ProngEnergyXView.push_back(EprongX);
386 ProngEnergyYView.push_back(EprongY);
392 double MaxEnergy1 = -999.0;
393 double MaxEnergy2 = -999.0;
395 if( MaxEnergy1<ProngEnergy[
i] ){
396 MaxEnergy1 = ProngEnergy[
i];
401 if(
i==Index1 )
continue;
402 if( MaxEnergy2<ProngEnergy[
i] ){
403 MaxEnergy2=ProngEnergy[
i];
407 double Eden = ProngEnergy[Index1]+ProngEnergy[Index2];
408 double Enum = ProngEnergy[Index1]-ProngEnergy[Index2];
409 if( Eden>0.0 ) Ebalance=Enum/Eden;
413 std::vector<int> Prong2DInXView;
414 std::vector<double> ProngEnergy2D;
415 int nprongs2D=SelectedProng2D.size();
416 for(
int ip=0;
ip<nprongs2D; ++
ip ){
419 int prongview=SelectedProng2D[
ip]->View();
420 Prong2DInXView.push_back(prongview);
425 const unsigned int nselectedprong2d_cells = current_prong->NCell();
426 for(
unsigned int icell=0; icell<nselectedprong2d_cells; ++icell){
427 const rb::RecoHit rhit = current_prong->RecoHit(icell);
430 Eprong += rhit.
GeV();
432 ProngEnergy2D.push_back(Eprong);
436 double max3D = -99.0;
438 if( max3D<ProngEnergy[
i] ){
439 max3D = ProngEnergy[
i];
445 for(
int i=0;
i<nprongs2D; ++
i ){
446 if( max2D<ProngEnergy2D[
i] ){
447 max2D = ProngEnergy2D[
i];
451 double ebalance2D = 1.0;
452 if( nprongs>0 && nprongs2D>0 ){
455 if( Prong2DInXView[index2D]==0 ){
456 Eden = ProngEnergy2D[index2D]+ProngEnergyXView[index3D];
457 Enum = ProngEnergy2D[index2D]-ProngEnergyXView[index3D];
460 Eden = ProngEnergy2D[index2D]+ProngEnergyYView[index3D];
461 Enum = ProngEnergy2D[index2D]-ProngEnergyYView[index3D];
463 if( Eden > 0.0 ) ebalance2D=
fabs(Enum)/Eden;
469 if (nCells > 0)
fInputVarsRVP[3] = (Float_t)longestTrack_ncells/nCells;
475 if (recoEnergy > 0.0)
fInputVarsRVP[6] = (Float_t)Eplane20/recoEnergy;
479 if (recoEnergy > 0.0)
fInputVarsRVP[9] = (Float_t)totalE_road2sigma/recoEnergy;
485 if (Esig_3 > 0.0)
fInputVarsRVP[14] = (Float_t)(recoEnergy - Esig_3)/Esig_3;
516 pidcol->push_back(pid);
525 evt.
put(std::move(pidcol));
526 evt.
put(std::move(assncol));
538 const unsigned int nvars = varnames_from_xml.size();
541 for(
unsigned int i=0;
i<
nvars; ++
i){
549 rvp.BookMVA(method_name, weight_file);
558 const bool get_full_name
564 cmd +=
"| grep \"<MethodSetup\"";
574 stream = popen(cmd.c_str(),
"r");
576 while ( fgets(buffer,
MAX_BUFFER, stream) != NULL ){
592 if(get_full_name)
return data;
597 std::vector<std::string>
names;
600 const boost::char_separator<char>
sep(
":");
601 boost::tokenizer< boost::char_separator<char> >
tokens(data, sep);
612 out_string = names[names.size() - 1];
623 cmd +=
"| grep \"Variable VarIndex=\" | awk \'{print $3}\'";
633 stream = popen(cmd.c_str(),
"r");
635 while ( fgets(buffer,
MAX_BUFFER, stream) != NULL ){
650 std::vector<std::string>
vars;
653 const boost::char_separator<char>
sep(
";");
654 boost::tokenizer< boost::char_separator<char> >
tokens(data, sep);
std::string fTrackLabel
Label for tracks.
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
void SetRecoE(float recE)
std::string getMethodNameFromXML(const std::string weights_file_name, const bool get_full_name=true) const
Get TMVA method name from XML.
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
fvar< T > fabs(const fvar< T > &x)
RecVarPID(fhicl::ParameterSet const &pset)
unsigned short Plane() const
void SetEFrac2SigRd(float eFrac2Sig)
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Vertical planes which measure X.
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
A collection of associated CellHits.
std::string fInstLabel3D
Instance label for 3d prongs.
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
void SetProngEBal3D(float pngEBal)
std::string fbdtRVP12MethodName
::xsd::cxx::tree::buffer< char > buffer
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
std::string fInstLabel2D
Instance label for 2d prongs.
const XML_Char const XML_Char * data
void SetNCells(unsigned int ncell)
ProductID put(std::unique_ptr< PROD > &&product)
virtual double TotalLength() const
Length (cm) of all the track segments.
bool replace_all(std::string &in, std::string const &from, std::string const &to)
Replace all occurrences of from in string with to.
TMVA::Reader * fbdtRVP12
Reader for RVP15.
void SetProngEBal2D(float pngEBal2D)
void SetProngs3D(int png3D)
static constexpr double L
void SetEFrac6PlWin(float eFrac6Plane)
bool fRunTMVA
Run the actual (memory-hungry) TMVA step?
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
void SetProngs2D(int png2D)
T get(std::string const &key) const
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
std::string fVertexLabel
Label for vertex.
void SetEIso3Sig(float eIso3Sig)
Vertex location in position and time.
void SetLongTr(float lTrack)
std::string fbdtRVPMethodName
TMVA::Reader * fbdtRVP
Reader for RVP15.
void produce(art::Event &evt)
void SetMipHits(float mHits)
void SetMipFrac(float mFrac)
unsigned int NYCell() const
Number of cells in the y-view.
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
std::string fRVP12WeightFile
RVP12 weight file name.
std::vector< std::string > getVarNamesFromXML(const std::string weights_file_name) const
Method to get Variable names from XML.
bool bookVariables(const std::string &weight_file, TMVA::Reader &rvp, std::string &method_name)
Book rvp reader.
#define MAX_BUFFER
for std::out
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
void SetEFrac2PlWin(float eFrac2Plane)
unsigned int NXCell() const
Number of cells in the x-view.
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
std::string fSliceLabel
Label for slices.
std::string fProngLabel
Label for prong.
void SetEPl20Frac(float ePlane20)
void SetLongTrFrac(float ltFrac)
bool IsNoise() const
Is the noise flag set?
std::string fRVPWeightFile
RVP weight file name.
rb::RecoHit MakeRecoHit(art::Ptr< rb::CellHit > const &chit, double const &w)
void SetRVP12(float rvp12)
std::vector< Float_t > fInputVarsRVP
Vector of the variables that are input for TMVA.