diff --git a/inc/VDispAnalyzer.h b/inc/VDispAnalyzer.h index 28504818..37381ac4 100644 --- a/inc/VDispAnalyzer.h +++ b/inc/VDispAnalyzer.h @@ -35,6 +35,8 @@ class VDispAnalyzer float fdistance_max; float floss_max; float fFui_min; + float fWidth_min; + float fFitstat_min; bool fDispErrorWeighting; float fDispErrorExponential; @@ -109,7 +111,8 @@ class VDispAnalyzer double* iR, double iEHeight, double iMCEnergy, double* img_fui, - float* img_pedvar ); + float* img_pedvar, + int* img_fitstat ); void calculateMeanDirection( float& xs, float& ys, vector< float >& x, vector< float >& y, @@ -134,7 +137,8 @@ class VDispAnalyzer double* img_fui, float* img_pedvar, double* pointing_dx, double* pointing_dy, - bool UseIntersectForHeadTail ); + bool UseIntersectForHeadTail, + int* img_fitstat ); vector< float > calculateExpectedDirectionError_or_Sign( unsigned int i_ntel, float iArrayElevation, float iArrayAzimuth, ULong64_t* iTelType, @@ -144,7 +148,7 @@ class VDispAnalyzer double* img_tgrad, double* img_loss, int* img_ntubes, double* img_weight, double xoff_4, double yoff_4, - double* img_fui, float* img_pedvar ); + double* img_fui, float* img_pedvar, int* img_fitstat ); float evaluate( float iWidth, float iLength, float iAsymm, float iDist, float iSize, float iPedvar, float itgrad, float iLoss, @@ -241,13 +245,15 @@ class VDispAnalyzer } void setQualityCuts( unsigned int iNImages_min = 0, float iAxesAngles_min = 0., float imaxdist = 1.e5, float imaxloss = 1., - float iminfui = 0. ) + float iminfui = 0., float iminwidth = -1., float iminfitstat = -10. ) { fAxesAngles_min = iAxesAngles_min; fNImages_min = iNImages_min; fdistance_max = imaxdist; floss_max = imaxloss; fFui_min = iminfui; + fWidth_min = iminwidth; + fFitstat_min = iminfitstat; } void setTelescopeTypeList( vector iTelescopeTypeList ); void setTelescopeFOV( vector< float > iTelFOV ) diff --git a/inc/VTableLookupRunParameter.h b/inc/VTableLookupRunParameter.h index 123e15a0..ded5e4ac 100644 --- a/inc/VTableLookupRunParameter.h +++ b/inc/VTableLookupRunParameter.h @@ -79,6 +79,8 @@ class VTableLookupRunParameter : public TNamed, public VGlobalRunParameter double fmaxloss; // note: same for all telescope types double fminfui ; // note: same for all telescope types double fminsize; // note: same for all telescope times + double fminwidth; // note: same for all telescope times + double fminfitstat; // note: same for all telescope times // seed for random selection of showers before table filling double fSelectRandom; int fSelectRandomSeed; @@ -110,6 +112,6 @@ class VTableLookupRunParameter : public TNamed, public VGlobalRunParameter void print( int iB = 0 ); void printHelp(); - ClassDef( VTableLookupRunParameter, 30 ); + ClassDef( VTableLookupRunParameter, 31 ); }; #endif diff --git a/src/VDispAnalyzer.cpp b/src/VDispAnalyzer.cpp index 83e24ee1..a3bab432 100644 --- a/src/VDispAnalyzer.cpp +++ b/src/VDispAnalyzer.cpp @@ -462,7 +462,8 @@ void VDispAnalyzer::calculateMeanDispDirection( unsigned int i_ntel, float* img_pedvar, double* pointing_dx, double* pointing_dy, - bool UseIntersectForHeadTail ) + bool UseIntersectForHeadTail, + int* img_fitstat ) { // reset values from previous event f_disp = -99.; @@ -477,7 +478,7 @@ void VDispAnalyzer::calculateMeanDispDirection( unsigned int i_ntel, || !img_asym || !img_tgrad || !img_loss || !img_ntubes || !img_weight || !img_fui - || !img_pedvar ) + || !img_pedvar || !img_fitstat ) { return; } @@ -501,7 +502,9 @@ void VDispAnalyzer::calculateMeanDispDirection( unsigned int i_ntel, if( img_size[i] > 0. && img_length[i] > 0. && img_ntubes[i] > 0 && sqrt( img_cen_x[i]*img_cen_x[i] + img_cen_y[i]*img_cen_y[i] ) < fdistance_max && img_loss[i] < floss_max - && img_fui[i] > fFui_min ) + && img_fui[i] > fFui_min + && img_width[i] > fWidth_min + && ( img_fitstat[i] < 1 || img_fitstat[i] >= fFitstat_min ) ) { disp = evaluate( ( float )img_width[i], ( float )img_length[i], ( float )img_asym[i], ( float )sqrt( img_cen_x[i] * img_cen_x[i] + img_cen_y[i] * img_cen_y[i] ), @@ -597,7 +600,8 @@ vector< float > VDispAnalyzer::calculateExpectedDirectionError_or_Sign( unsigned double xoff_4, double yoff_4, double* img_fui, - float* img_pedvar ) + float* img_pedvar, + int* img_fitstat ) { vector< float > i_disp( i_ntel, -9999. ); // make sure that all data arrays exist @@ -607,7 +611,7 @@ vector< float > VDispAnalyzer::calculateExpectedDirectionError_or_Sign( unsigned || !img_asym || !img_tgrad || !img_loss || !img_ntubes || !img_weight || !img_fui - || !img_pedvar ) + || !img_pedvar || !img_fitstat ) { return i_disp; } @@ -621,7 +625,9 @@ vector< float > VDispAnalyzer::calculateExpectedDirectionError_or_Sign( unsigned if( img_size[i] > 0. && img_length[i] > 0. && sqrt( img_cen_x[i]*img_cen_x[i] + img_cen_y[i]*img_cen_y[i] ) < fdistance_max && img_loss[i] < floss_max - && img_fui[i] > fFui_min ) + && img_fui[i] > fFui_min + && img_width[i] > fWidth_min + && ( img_fitstat[i] < 1 || img_fitstat[i] >= fFitstat_min ) ) { i_disp[i] = evaluate( ( float )img_width[i], ( float )img_length[i], ( float )img_asym[i], ( float )sqrt( img_cen_x[i] * img_cen_x[i] + img_cen_y[i] * img_cen_y[i] ), @@ -706,7 +712,8 @@ void VDispAnalyzer::calculateEnergies( unsigned int i_ntel, double iEHeight, double iMCEnergy, double* img_fui, - float* img_pedvar ) + float* img_pedvar, + int* img_fitstat ) { fdisp_energy = -99.; fdisp_energy_chi = 1.e-10; @@ -723,7 +730,7 @@ void VDispAnalyzer::calculateEnergies( unsigned int i_ntel, || !img_asym || !img_tgrad || !img_loss || !img_ntubes || !img_weight || !iRcore - || !img_fui ) + || !img_fui || !img_fitstat ) { return; } @@ -738,7 +745,9 @@ void VDispAnalyzer::calculateEnergies( unsigned int i_ntel, if( img_size[i] > 0. && iRcore[i] > 0. && iArrayElevation > 0. && sqrt( img_cen_x[i]*img_cen_x[i] + img_cen_y[i]*img_cen_y[i] ) < fdistance_max && img_loss[i] < floss_max - && img_fui[i] > fFui_min ) + && img_fui[i] > fFui_min + && img_width[i] > fWidth_min + && ( img_fitstat[i] < 1 || img_fitstat[i] >= fFitstat_min ) ) { fdisp_energy_T[i] = fTMVADispAnalyzer->evaluate( ( float )img_width[i], ( float )img_length[i], diff --git a/src/VTableLookupDataHandler.cpp b/src/VTableLookupDataHandler.cpp index 61d0b2f9..d81a7eda 100644 --- a/src/VTableLookupDataHandler.cpp +++ b/src/VTableLookupDataHandler.cpp @@ -561,6 +561,7 @@ int VTableLookupDataHandler::fillNextEvent( bool bShort ) } } fCurrentNoiseLevel[i] = ftpars[i]->meanPedvar_Image; + fFitstat[i] = ftpars[i]->Fitstat; if( !bShort ) { fmeanPedvar_ImageT[i] = ftpars[i]->meanPedvar_Image; @@ -577,7 +578,6 @@ int VTableLookupDataHandler::fillNextEvent( bool bShort ) fmaxindex2[i] = ftpars[i]->index_of_max[1]; fmaxindex3[i] = ftpars[i]->index_of_max[2]; ftchisq_x[i] = ftpars[i]->tchisq_x; - fFitstat[i] = ftpars[i]->Fitstat; } } else @@ -616,7 +616,9 @@ int VTableLookupDataHandler::fillNextEvent( bool bShort ) fDispAnalyzerEnergy->setQualityCuts( fSSR_NImages_min, fSSR_AxesAngles_min, fTLRunParameter->fmaxdist, fTLRunParameter->fmaxloss, - fTLRunParameter->fminfui ); + fTLRunParameter->fminfui, + fTLRunParameter->fminwidth, + fTLRunParameter->fminfitstat ); fDispAnalyzerEnergy->calculateEnergies( getNTel(), fArrayPointing_Elevation, fArrayPointing_Azimuth, @@ -632,7 +634,7 @@ int VTableLookupDataHandler::fillNextEvent( bool bShort ) getDistanceToCoreTel(), fEmissionHeightMean, fMCEnergy, - ffui, fmeanPedvar_ImageT ); + ffui, fmeanPedvar_ImageT, fFitstat ); // fill results setEnergy( fDispAnalyzerEnergy->getEnergy(), false ); @@ -718,7 +720,7 @@ void VTableLookupDataHandler::doStereoReconstruction() floss, fntubes, getWeight(), fXoff_intersect, fYoff_intersect, - ffui, fmeanPedvar_ImageT ); + ffui, fmeanPedvar_ImageT, fFitstat ); } //////////////////////////////////////////////////////////////////// // estimate disp head/tail sign @@ -738,7 +740,7 @@ void VTableLookupDataHandler::doStereoReconstruction() floss, fntubes, getWeight(), fXoff_intersect, fYoff_intersect, - ffui, fmeanPedvar_ImageT ); + ffui, fmeanPedvar_ImageT, fFitstat ); } // use weighting calculated from disp error @@ -747,7 +749,9 @@ void VTableLookupDataHandler::doStereoReconstruction() fDispAnalyzerDirection->setQualityCuts( fSSR_NImages_min, fSSR_AxesAngles_min, fTLRunParameter->fmaxdist, fTLRunParameter->fmaxloss, - fTLRunParameter->fmaxloss ); + fTLRunParameter->fminfui, + fTLRunParameter->fminwidth, + fTLRunParameter->fminfitstat ); fDispAnalyzerDirection->setTelescopeFOV( fTelFOV ); fDispAnalyzerDirection->calculateMeanDispDirection( getNTel(), @@ -764,7 +768,7 @@ void VTableLookupDataHandler::doStereoReconstruction() iDispError, iDispSign, ffui, fmeanPedvar_ImageT, fpointing_dx, fpointing_dy, - fTLRunParameter->fDisp_UseIntersectForHeadTail ); + fTLRunParameter->fDisp_UseIntersectForHeadTail, fFitstat ); // reconstructed direction by disp method: fXoff = fDispAnalyzerDirection->getXcoordinate_disp(); fYoff = fDispAnalyzerDirection->getYcoordinate_disp(); @@ -1487,6 +1491,8 @@ bool VTableLookupDataHandler::setOutputFile( string iOutput, string iOption, str fOTree->Branch( "ntubes", fntubes, iTT ); sprintf( iTT, "nsat[%d]/s", fNTel ); fOTree->Branch( "nsat", fnsat, iTT ); + sprintf( iTT, "fui[%d]/D", fNTel ); + fOTree->Branch( "fui", ffui, iTT ); sprintf( iTT, "nlowgain[%d]/s", fNTel ); fOTree->Branch( "nlowgain", fnlowgain, iTT ); sprintf( iTT, "alpha[%d]/D", fNTel ); @@ -1507,9 +1513,9 @@ bool VTableLookupDataHandler::setOutputFile( string iOutput, string iOption, str fOTree->Branch( "tgrad_x", ftgrad_x, iTT ); sprintf( iTT, "tchisq_x[%d]/D", fNTel ); fOTree->Branch( "tchisq_x", ftchisq_x, iTT ); - sprintf( iTT, "Fitstat[%d]/I", fNTel ); - fOTree->Branch( "Fitstat", fFitstat, iTT ); } + sprintf( iTT, "Fitstat[%d]/I", fNTel ); + fOTree->Branch( "Fitstat", fFitstat, iTT ); fOTree->Branch( "DispNImages", &fnxyoff, "DispNImages/i" ); fOTree->Branch( "DispXoff_T", fXoff_T, "DispXoff_T[NImages]/F" ); fOTree->Branch( "DispYoff_T", fYoff_T, "DispYoff_T[NImages]/F" ); diff --git a/src/VTableLookupRunParameter.cpp b/src/VTableLookupRunParameter.cpp index 9809779b..6c05af1d 100644 --- a/src/VTableLookupRunParameter.cpp +++ b/src/VTableLookupRunParameter.cpp @@ -38,6 +38,8 @@ VTableLookupRunParameter::VTableLookupRunParameter() fmaxdist = 50000.; fmaxloss = 1.; fminfui = 0.; + fminwidth = -1.; + fminfitstat = -10.; fSelectRandom = -1.; fSelectRandomSeed = 17; fRerunStereoReconstruction = false; @@ -294,6 +296,14 @@ bool VTableLookupRunParameter::fillParameters( int argc, char* argv[] ) { fminfui = atof( iTemp.substr( iTemp.rfind( "=" ) + 1, iTemp.size() ).c_str() ); } + else if( iTemp.find( "-minwidth" ) < iTemp.size() ) + { + fminwidth = atof( iTemp.substr( iTemp.rfind( "=" ) + 1, iTemp.size() ).c_str() ); + } + else if( iTemp.find( "-minfitstat" ) < iTemp.size() ) + { + fminfitstat = atof( iTemp.substr( iTemp.rfind( "=" ) + 1, iTemp.size() ).c_str() ); + } else if( iTemp.find( "-maxdistancetocameracenter" ) < iTemp.size() ) { fMC_distance_to_cameracenter_max = atof( iTemp.substr( iTemp.rfind( "=" ) + 1, iTemp.size() ).c_str() ); @@ -564,7 +574,15 @@ void VTableLookupRunParameter::print( int iP ) } if( fminfui > 0. ) { - cout << "\t BDT TMVA stereo reconstruction fui cut < " << fminfui << endl; + cout << "\t BDT TMVA stereo reconstruction fui cut > " << fminfui << endl; + } + if( fminwidth > 0. ) + { + cout << "\t BDT TMVA stereo reconstruction min width cut < " << fminwidth << endl; + } + if( fminfitstat > 0. ) + { + cout << "\t BDT TMVA stereo reconstruction min fitstat cut < " << fminfitstat << endl; } cout << "\t Head/tail uncertainty: "; if( fDisp_UseIntersectForHeadTail )