RootCore Packages: /cvmfs/atlas.cern.ch/repo/sw/ASG/AnalysisBase/2.4.31/DiTauMassTools/DiTauMassTools/MissingMassCalculator.h Source File

00001 /*
00002 >>>>>> MissingMassCalculator
00003 >>>>>> Author: Aliaksandr Pranko (appranko@lbl.gov)
00004 >>>>>> Code developers: David Rousseau (rousseau@lal.in2p3.fr), Dimitris Varouchas (Dimitris.Varouchas@cern.ch)
00005 MissingMassCalculator is designed to reconstruct mass in
00006 events where two particles decay into states with missing ET.
00007 */
00008 #ifndef MissingMassCalculator_h
00009 #define MissingMassCalculator_h
00010 
00011 
00012 
00013 #if !defined (__CINT__) || defined (__MAKECINT__)
00014 
00015 #include <TH1.h>
00016 #include <TF1.h>
00017 #include <TMath.h>
00018 #include <TLorentzVector.h>
00019 #include <TVector3.h>
00020 #include <vector>
00021 #include <TObject.h>
00022 #include <TDirectory.h>
00023 
00024 #include <memory>
00025 
00026 #include <TRandom2.h>
00027 
00028 class TH1F;
00029 
00030 
00031 namespace MMCCalibrationSet 
00032 {
00033   enum e { MMC2011=0, MMC2012, MMC2015, MMC2015HIGHMASS, UPGRADE, LFVMMC2012, MMC2016MC15C, MAXMMCCALIBRATIONSET };
00034   const std::string name[MAXMMCCALIBRATIONSET]={ "MMC2011", "MMC2012", "MMC2015", "MMC2015HIGHMASS", "UPGRADE", "LFVMMC2012", "MMC2016MC15C"};
00035 }
00036 
00037 namespace MMCFitMethod 
00038 {
00039   enum e { MAXW=0, MLM, MLNU3P,MAX};
00040   const std::string name[MAX]={ "MAXW=MaximumWeight", "MLM=MostLikelyMass", "MLNU3P=MostLikelyNeUtrino3Momentum"};
00041   const std::string shortName[MAX]={ "MAXW", "MLM", "MLNU3P"};
00042 }
00043 
00044 namespace WalkStrategy 
00045 {
00046   enum e { GRID=0, RANDOM,RANDOMNONUNIF, MARKOVCHAIN, MAXWALKSTRATEGY };
00047 }
00048 
00049 namespace MaxHistStrategy 
00050 {
00051   enum e { MAXBIN=0,MAXBINWINDOW, SLIDINGWINDOW, FIT,MAXMAXHISTSTRATEGY };
00052 }
00053 
00054 
00055 namespace HistInfo 
00056 {
00057   enum e { PROB=0, INTEGRAL, CHI2, DISCRI, TANTHETA, TANTHETAW, FITLENGTH, RMS, RMSVSDISCRI, MEANBIN, MAXHISTINFO };
00058 }
00059 
00060 
00061 
00062 #endif
00063 
00064 /* class MissingMassCalculator: public TObject { */
00065 class MissingMassCalculator {
00066 
00067  public:
00068   //DR
00069   /* #if !defined(__CINT__) | defined(__MAKECINT__) */
00070   /* // Define the class for the cint dictionary */
00071   /* ClassDef(MissingMassCalculator,1); */
00072   /* #endif */
00073     
00074 
00075   // keep the whole list of points
00076   struct AllPointsOutputInfo {
00077     std::vector<double> prob;
00078     std::vector<double> mtautau;
00079     std::vector<double> Phi1;
00080     std::vector<double> Phi2;
00081     std::vector<double> MEtP;
00082     std::vector<double> MEtL;
00083     std::vector<double> MEtX;
00084     std::vector<double> MEtY;
00085     std::vector<double> Mnu1;
00086     std::vector<double> Mnu2;
00087     std::vector<double> eTau1;
00088     std::vector<double> eTau2;
00089     std::vector<int> pCode; // -10 : outside parameter space boundary. 0 no solution. 10 some solutions. +1 if accepted
00090     std::vector<int> index;
00091   };
00092 
00093   
00094       
00095 
00096 private:
00097   AllPointsOutputInfo m_aPOI;
00098   int m_apoiIndex;
00099 
00100   //---------------- structures
00101   struct DitauStuff {
00102     double Mditau_best; // best fitted M(ditau)
00103     double Sign_best;   // best significance of M(ditau) fit
00104     TLorentzVector nutau1;  // fitted 4-vec for neutrino from tau-1
00105     TLorentzVector nutau2;  // fitted 4-vec for neutrino from tau-2
00106     TLorentzVector vistau1;  // fitted 4-vec for visible tau-1
00107     TLorentzVector vistau2;  // fitted 4-vec for visible tau-2
00108     TLorentzVector taufit1;  // fitted 4-vec for full tau-1
00109     TLorentzVector taufit2;  // fitted 4-vec for full tau-2
00110     double RMSoverMPV;
00111   };
00112 
00113   struct InputInfoStuff {
00114     TVector2 MetVec;
00115     TVector2 detMetVec; // need to keep separate instead Metvec recalibrated or replaced by -HT
00116     TLorentzVector vistau1;
00117     TLorentzVector vistau2;
00118     int type_visTau1; // 0: l, 1:1-prong, 3:3-prong
00119     int type_visTau2;  // 0: l, 1:1-prong, 3:3-prong
00120     int Nprong_tau1;
00121     int Nprong_tau2;
00122     int dataType;
00123     double METcovphi;
00124     double METsigmaP;
00125     double METsigmaL;
00126     double SumEt;
00127     double sigmaEtau1;
00128     double sigmaEtau2;
00129     std::vector<TLorentzVector> jet4vecs;
00130     int Njet25;
00131     double DelPhiTT;
00132     bool allowUseHT;
00133     bool UseHT;
00134     double MHtSigma1; // sigma of 1st Gaussian in missing Ht resolution
00135     double MHtSigma2; // sigma of 2nd Gaussian in missing Ht resolution
00136     double MHtGaussFr; // relative fraction of 2nd Gaussian
00137     double HtOffset; // HT offset
00138 
00139   };
00140 
00141   struct OutputInfoStuff {
00142     int FitStatus;
00143     double FitSignificance[MMCFitMethod::MAX];
00144     double FittedMass[MMCFitMethod::MAX];
00145     TLorentzVector nuvec1[MMCFitMethod::MAX];
00146     TLorentzVector objvec1[MMCFitMethod::MAX];
00147     TLorentzVector nuvec2[MMCFitMethod::MAX];
00148     TLorentzVector objvec2[MMCFitMethod::MAX];
00149     TLorentzVector totalvec[MMCFitMethod::MAX];
00150     TVector2 FittedMetVec[MMCFitMethod::MAX];
00151     double RMS2MPV;
00152     std::shared_ptr<TH1F> hMfit_all;
00153     std::shared_ptr<TH1F> hMfit_allNoWeight;
00154     int NTrials;
00155     int NSuccesses;
00156     int NSolutions;
00157     double SumW;
00158     double AveSolRMS;
00159     bool UseHT;
00160   };
00161 
00162   MMCCalibrationSet::e m_MMCCalibrationSet;
00163   //  WalkStrategy::e walkStrategy;
00164   TRandom2 * randomGen;
00165   
00166   //SpeedUp static array for efficient access
00167   static double fit_param[2][3][6][5];
00168   static double ter_sigma_par[2][10][3];
00169   //cache quantities for efficient NuPSolution calculation
00170 
00171 
00172   int nsolmax,nsolfinalmax;
00173   int m_NiterRandomLocal;
00174   int m_NsucStop;
00175   int m_rmsStop;
00176   double m_meanbinStop;
00177   
00178 
00179 
00180   
00181   // these are temporary vectors. Dclared globally to avoid construction/destruction
00182   std::vector<TLorentzVector> nuvecsol1;
00183   std::vector<TLorentzVector> nuvecsol2;
00184 
00185   std::vector<TLorentzVector> tauvecsol1;
00186   std::vector<TLorentzVector> tauvecsol2;
00187   std::vector<double> tauvecprob1;
00188   std::vector<double> tauvecprob2;
00189 
00190   std::vector<TLorentzVector> nuvec1_tmp;
00191   std::vector<TLorentzVector> nuvec2_tmp;
00192 
00193   TVector2 m_metVec;
00194   TLorentzVector tautau_tmp;
00195   TLorentzVector tlv_tmp;
00196 
00197   
00198   //  std::vector<TLorentzVector> tauvecsol1;
00199   // std::vector<TLorentzVector> tauvecsol2;
00200   bool debugThisIteration;
00201   
00202   int nCallprobCalculatorV9fast;
00203   
00204 
00205   int iter1,iter2,iter3,iter4,iter5,iang1low,iang1high,iang2low,iang2high;
00206   int iterTheta3d;
00207   
00208   double prob_tmp;
00209   
00210   double totalProbSum;
00211   double mtautauSum;
00212 
00213   int m_seed;
00214   // data member for the spaceWalker approach
00215 
00216   int m_iter0;
00217   int m_iterNuPV3;
00218   int m_testptn1;
00219   int m_testptn2;
00220   int m_testdiscri1;
00221   int m_testdiscri2;
00222   int m_iterNsuc;
00223   bool m_switch1;
00224   bool m_switch2;
00225 
00226   bool m_meanbinToBeEvaluated;
00227   
00228   int m_markovCountDuplicate;
00229   int m_markovNFullScan;
00230   int m_markovNRejectNoSol;
00231   int m_markovNRejectMetropolis;
00232   int m_markovNAccept;
00233 
00234   double m_PrintmMaxError;
00235   double m_PrintmMeanError;
00236   double m_PrintmInvWidth2Error;
00237 
00238   double m_ProposalTryMEt;
00239   double m_ProposalTryPhi;
00240   double m_ProposalTryMnu;
00241   double m_ProposalTryEtau;
00242 
00243 
00244   double m_mTau,m_mTau2;
00245   double m_MEtL,m_MEtP,m_Phi1,m_Phi2,m_Mnu1,m_Mnu2;
00246   double m_eTau1, m_eTau2;
00247   double m_eTau10, m_eTau20;
00248   double m_MEtL0,m_MEtP0,m_Phi10,m_Phi20,m_Mnu10,m_Mnu20;
00249   double m_MEtLMin,m_MEtPMin,m_Phi1Min,m_Phi2Min,m_Mnu1Min,m_Mnu2Min;
00250   double m_MEtLMax,m_MEtPMax,m_Phi1Max,m_Phi2Max,m_Mnu1Max,m_Mnu2Max;
00251   double m_MEtLStep,m_MEtPStep,m_Phi1Step,m_Phi2Step,m_Mnu1Step,m_Mnu2Step;
00252   double m_MEtLRange,m_MEtPRange,m_Phi1Range,m_Phi2Range,m_Mnu1Range,m_Mnu2Range;
00253   double m_MEtProposal,m_PhiProposal,m_MnuProposal;
00254   double m_eTau1Proposal,m_eTau2Proposal;
00255 
00256   double m_MEtX, m_MEtY,m_MEtT;
00257   double m_eTau1Min,m_eTau1Max,m_eTau1Range;
00258   double m_eTau2Min,m_eTau2Max,m_eTau2Range;
00259   bool fullParamSpaceScan;
00260   bool m_Mnu1Exclude;
00261   int m_nsolOld;
00262   std::vector<double>   m_probFinalSolOldVec;
00263   std::vector<double>   m_mtautauFinalSolOldVec;
00264   std::vector<TLorentzVector>  m_nu1FinalSolOldVec;
00265   std::vector<TLorentzVector>  m_nu2FinalSolOldVec;
00266 
00267   int m_nsol;
00268   std::vector<double>   m_probFinalSolVec;
00269   std::vector<double>   m_mtautauFinalSolVec;
00270   std::vector<TLorentzVector>  m_nu1FinalSolVec;
00271   std::vector<TLorentzVector>  m_nu2FinalSolVec;
00272 
00273 
00274   double m_Mnu1ExcludeMin,m_Mnu1ExcludeMax,m_Mnu1ExcludeRange;
00275   double m_Mnu1XMin,m_Mnu1XMax, m_Mnu1XRange;
00276   double m_walkWeight;
00277   double m_cosPhi1, m_cosPhi2, m_sinPhi1, m_sinPhi2;
00278   
00279   bool m_scanMnu1,m_scanMnu2;
00280   bool m_scanEtau1,m_scanEtau2;
00281 
00282   TLorentzVector m_tauVec1,m_tauVec2;
00283   double m_tauVec1Phi, m_tauVec2Phi;
00284   double m_tauVec1M, m_tauVec2M;
00285   double m_tauVec1Px, m_tauVec1Py, m_tauVec1Pz;
00286   double m_tauVec2Px, m_tauVec2Py, m_tauVec2Pz;
00287   double m_tauVec1P, m_tauVec2P;
00288   double m_tauVec1E;
00289   double m_tauVec2E;
00290   double m_m2Nu1;
00291   double m_m2Nu2;
00292   double m_ET2v1;
00293   double m_ET2v2;
00294   double m_E2v1;
00295   double m_E2v2;
00296   double m_Ev2;
00297   double m_Ev1;
00298   double m_metCovPhiCos,m_metCovPhiSin;
00299   double m_Mvis,m_Meff;
00300   double m_DetMEt;
00301   double m_inputMEtX,m_inputMEtY,m_inputMEtT;
00302   double m_HtOffset;
00303   bool reRunWithBestMET;
00304   
00305   //--- define histograms for histogram method
00306   //--- upper limits need to be revisied in the future!!! It may be not enough for some analyses
00307   std::shared_ptr<TH1F> fMfit_all;
00308   std::shared_ptr<TH1F> fMfit_allNoWeight;
00309   
00310   std::shared_ptr<TH1F> fPXfit1;
00311   std::shared_ptr<TH1F> fPYfit1;
00312   std::shared_ptr<TH1F> fPZfit1;
00313   std::shared_ptr<TH1F> fPXfit2;
00314   std::shared_ptr<TH1F> fPYfit2;
00315   std::shared_ptr<TH1F> fPZfit2;
00316 
00317   TF1 *m_fFitting;
00318 
00319   TH1F* fPhi1;
00320   TH1F* fPhi2;
00321   TH1F* fMnu1;
00322   TH1F* fMnu2;
00323   TH1F* fMetx;
00324   TH1F* fMety;
00325   TH1F* fTheta3D;
00326   TH1F* fTauProb;
00327 
00328   // for intermediate calc
00329   TLorentzVector TLVdummy;
00330 
00331   //---------------- protected variables
00332   DitauStuff fDitauStuffFit; // results based on fit method
00333   DitauStuff fDitauStuffHisto; // results based on histo method
00334   InputInfoStuff rawInput; // raw input 
00335   InputInfoStuff preparedInput; // corrected input (DR should maybe not be a datamember)
00336   OutputInfoStuff OutputInfo; // output info
00337 
00338   int fUseVerbose; // code to turn ON printouts for debugging
00339   bool fSpeedStudy; // code to turn ON speed study
00340   int AlgorithmVersion; // version of the algorithm
00341   int SearchMode; // search Mode: 0=di-tau, 1=WW, 2=W->tau+nu
00342   int fApplyMassScale; // switch to apply mass scale correction
00343   int fUseTailCleanup; // switch to apply tail clean-up
00344   int fUseTauProbability; // switch to apply TauProbability
00345   int fUseMnuProbability; // switch to apply MnuProbability
00346   int fUseDefaults; // switch to control defaults: 1== use defaults, 0== don't use defaults (useful for studies) 
00347   int fUseEfficiencyRecovery; // switch to turn ON/OFF re-fit in order to recover efficiency
00348 
00349   int Niter_fit1; // number of iterations for dR-dPhi scan 
00350   int Niter_fit2; // number of iterations for MET-scan 
00351   int Niter_fit3; // number of iterations for Mnu-scan 
00352   int NiterRandom; // number of random iterations (for lh, multiply or divide by 10 for ll and hh)
00353   int NsucStop;
00354   int RMSStop;
00355   int RndmSeedAltering; // reset seed (not necessary by default)
00356 
00357   int InputReorder; // flag for input re-order
00358   int LFVmode; // flag to determine which LFV decay to be reconstructed: 0=H->e+tau(mu) or 1=H->mu+tau(e) 
00359 
00360   int fJERsyst; // switch for JER systematics
00361   int METresSyst; // switch to turn on/off MET resolution systematics
00362   double dTheta3d_binMin; // minimal step size for dTheta3D
00363   double dTheta3d_binMax; // maximum step size for dTheta3D
00364   double dRmax_tau; // maximum dR(nu-visTau)
00365   double Nsigma_METscan,Nsigma_METscan2,Nsigma_METscan_ll,Nsigma_METscan_lh,Nsigma_METscan_hh; // number of sigmas for MET-scan
00366 
00367   double beamEnergy; // beam energy (Tevatron=980, LHC-1=3500.0) 
00368   int METScanScheme; // MET-scan scheme: 0- use JER; 1- use simple sumEt & missingHt for Njet=0 events in (lep-had) 
00369   double MnuScanRange; // range of M(nunu) scan; M(nunu) range can be affected by selection cuts
00370   
00371   bool fUseDphiLL; //for leplep
00372 
00373   //---------------- protected functions
00374   void ClearDitauStuff(DitauStuff &fStuff);
00375   void ClearInputStuff();
00376   void ClearOutputStuff();
00377   void DoOutputInfo();
00378   void DoMetResolution(InputInfoStuff &fStuff);
00379   void PrintInputInfo(const InputInfoStuff & fStuff);
00380   void PrintOtherInput();
00381   void PrintResults();
00382   void FinalizeInputStuff();//const InputInfoStuff &rawInput, InputInfoStuff &preparedInput);
00383   int NuPsolution(TVector2 met_vec, double theta1, double phi1, 
00384           double theta2, double phi2, double &P1, double &P2); // keep this version for simple tests
00385 
00386 
00387   inline int NuPsolutionV3(const double & mNu1, const double & mNu2, const double & phi1, const double & phi2, 
00388                int & nsol1, int & nsol2);
00389 
00390   inline int NuPsolutionLFV(const TVector2 & met_vec, const TLorentzVector & tau, 
00391                 const double & m_nu, std::vector<TLorentzVector> &nu_vec);
00392 
00393   
00394   inline double tauSigmaE(const TLorentzVector & tauVec, const int tauType );
00395   inline double dTheta3Dparam(const int & parInd, const int & tau_type, const double & P_tau,const double *par); 
00396 
00397   //DR public:  
00398   double dTheta3d_probabilityFast(const int & tau_type,const double & dTheta3d,const double & P_tau);
00399  protected:
00400   inline double myDelThetaHadFunc(double *x, double *par);
00401   inline double myDelThetaLepFunc(double *x, double *par);
00402 
00403   inline double MetProbability(const double & met1,const  double & met2,const  double & MetSigma1, const double & MetSigma2);
00404 
00405   inline double MHtProbability(const double & d_mhtX, const double & d_mhtY, const double & mht, 
00406                    const double & trueMetGuess, const double & mht_offset);
00407   inline double MHtProbabilityHH(const double & d_mhtX, const double & d_mhtY, const double & mht, 
00408                  const double & trueMetGuess, const double & mht_offset);
00409 
00410   inline double mEtAndTauProbability();
00411 
00412   inline double TERProbabilitySimple(int tau_type, const TLorentzVector & vec, 
00413                      const double & Pt_scan);
00414   inline double TERSimpleSigma(int tau_type, const TLorentzVector & vec);
00415 
00416   inline void swapTLVPointers (TLorentzVector * & p1,TLorentzVector * & p2);
00417   inline double Angle(const TLorentzVector & vec1, const TLorentzVector & vec2);
00418   inline double AngularProbability(TLorentzVector nu_vec, TLorentzVector vis_vec, int decayType);
00419   inline int CheckSolutions(TLorentzVector nu_vec, TLorentzVector vis_vec, int decayType);
00420   inline double MnuProbability(double mnu, double binsize);
00421   inline double MnuProbability(double mnu);
00422   inline double ResonanceProbability(TLorentzVector vec1, TLorentzVector vec2);
00423   inline int TailCleanUp(const int & type1, const TLorentzVector & vis1, const TLorentzVector & nu1, 
00424              const int & type2, const TLorentzVector & vis2, const TLorentzVector & nu2, 
00425              const double & mmc_mass, const double & vis_mass, const double & eff_mass, const double & dphiTT);
00426 
00427   inline double TauProbability(const int & type1, const TLorentzVector & vis1, const TLorentzVector & nu1, 
00428                    const int & type2, const TLorentzVector & vis2, const TLorentzVector & nu2);
00429   inline double TauProbability(const int & type1, const TLorentzVector & vis1, const TLorentzVector & nu1, 
00430                    const int & type2, const TLorentzVector & vis2, const TLorentzVector & nu2, const double & detmet);
00431   inline double TauProbabilityMatrix(const int & type1, const TLorentzVector & vis1, const TLorentzVector & nu1, 
00432                      const int & type2, const TLorentzVector & vis2, const TLorentzVector & nu2);
00433   inline double TauProbabilityLFV(const int & type1, const TLorentzVector & vis1, const TLorentzVector & nu1);
00434   inline double MaxDelPhi(int tau_type, double Pvis);
00435 
00436   inline int refineSolutions (          const double & M_nu1, const double & M_nu2,
00437                         const int nsol1, const int nsol2,
00438                   const double & Mvis, const double & Meff);
00439 
00440   inline double mT(const TLorentzVector & vec,const TVector2 & met_vec);
00441   
00442   
00443   inline void handleSolutions();
00444 
00445   inline double MassScale(int method, double mass, const int & tau_type1, const int & tau_type2);
00446 
00447   
00448   // factor out parameter space walking and probablity computing
00449   inline  int DitauMassCalculatorV9walk();
00450   
00451 
00452   // Calculates mass of lep+tau system in LFV X->lep+tau decays
00453 
00454   // It is based on DitauMassCalculatorV9, not optimized for speed yet, simple phase-space scan   
00455   inline int DitauMassCalculatorV9lfv();
00456 
00457 
00458   
00459   // only compute probability
00460   inline int probCalculatorV9fast(
00461                   const double & phi1, const double & phi2, 
00462                   const double & M_nu1, const double & M_nu2);
00463   
00464   // initialize the walker
00465   inline  void SpaceWalkerInit();
00466   
00467   //walk the walker
00468   inline bool SpaceWalkerWalk();
00469   
00470   //update tau kinematics if tau energy scan 
00471   inline void updateTauKine();
00472 
00473 
00474   inline bool precomputeCache();
00475   
00476 
00477   inline bool updateDouble  (const double in, double & out) ;
00478   inline bool checkMEtInRange  () ;
00479   inline bool checkAllParamInRange  () ;
00480   inline double fixPhiRange (const double & phi);
00481 // fast (but approximate) calculation of sin and cos
00482   inline void fastSinCos (const double & phi, double & sinPhi, double & cosPhi);
00483   
00484 
00485   //----------------------------------------------
00486   //
00487   // >>>>>>>>>>>>>   Public methods <<<<<<<<<<<< 
00488   //
00489   //______________________________________________
00490 
00491 public:
00492 
00493   ~MissingMassCalculator() ;
00494 
00495   MissingMassCalculator() ;
00496 
00497 
00498   int RunMissingMassCalculator();
00499   
00500   //-------- Set Input Parameters
00501   void SetCalibrationSet( MMCCalibrationSet::e aset) { m_MMCCalibrationSet=aset; }
00502   void SetSearchMode(int val) { SearchMode=val; }
00503   void SetUseVerbose(int val) { fUseVerbose=val; }
00504   void SetSpeedStudy(int val) { fSpeedStudy=val; }
00505   void SetAlgorithmVersion(int val) { AlgorithmVersion=val; }
00506   void SetUseTauProbability(int val) { fUseTauProbability=val; }
00507   void SetUseMnuProbability(int val) { fUseMnuProbability=val; }
00508   void SetUseDefaults(int val) { fUseDefaults=val; }
00509   void SetUseEfficiencyRecovery(int val) { fUseEfficiencyRecovery=val; }
00510 
00511   void SetNiterFit1(int val) { Niter_fit1=val; } // number of iterations per loop in dPhi loop
00512   void SetNiterFit2(int val) { Niter_fit2=val; } // number of iterations per loop in MET loop
00513   void SetNiterFit3(int val) { Niter_fit3=val; } // number of iterations per loop in Mnu loop
00514   void SetNiterRandom(int val) { NiterRandom=val; } // number of random iterations
00515   void SetNsucStop(int val) { NsucStop=val; } // Arrest criteria for Nsuccesses
00516   void SetRMSStop(int val) { RMSStop=val;}
00517   void SetMeanbinStop(double val) {m_meanbinStop=val;}
00518   void SetRndmSeedAltering(int val) { RndmSeedAltering=val; } // number of iterations per loop in Mnu loop
00519   void SetMaxDRtau(double val) { dRmax_tau=val; } // max value of dR
00520   void SetNsigmaMETscan(double val) { Nsigma_METscan=val; } // number of sigma's for MET-scan
00521   void SetNsigmaMETscan_ll(double val) { Nsigma_METscan_ll=val; } // number of sigma's for MET-scan in ll events
00522   void SetNsigmaMETscan_lh(double val) { Nsigma_METscan_lh=val; } // number of sigma's for MET-scan in lh events
00523   void SetNsigmaMETscan_hh(double val) { Nsigma_METscan_hh=val; } // number of sigma's for MET-scan in hh events
00524   void SetBeamEnergy(double val) { beamEnergy=val; } // beam energy
00525   void SetdTheta3d_binMax(double val) { dTheta3d_binMax=val; } // maximum step size for dTheta3D
00526   void SetdTheta3d_binMin(double val) { dTheta3d_binMin=val; } // minimal step size for dTheta3D
00527 
00528   void SetMetVec(const TVector2 & vec);
00529   void SetVisTauVec(int i, const TLorentzVector & vec);
00530   void SetVisTauType(int i, int tautype);
00531   void SetNprong(int i, int nprong);
00532   void SetSumEt(double sumEt);
00533   void SetIsData(int val);
00534   void SetMetCovariance(double varX, double varY, double varXY); // MET covariance matrix
00535   void SetMetScanParams(double phi, double sigmaP, double sigmaL); // MET ellipse (deduced from MetCovariance)
00536   void SetMetScanParamsUE(double sumEt, double phi_scan=0.0, int data_code=0);
00537   void SetMetScanParamsJets(std::vector<TLorentzVector> jets);
00538   void SetJERsyst(int val) { fJERsyst=val; }
00539   void SetMETresSyst(int val) { METresSyst=val; } // MET resolution systematics: +/-1 sigma
00540   void SetApplyMassScale(int val) { fApplyMassScale=val; } 
00541   void SetUseTailCleanup(int val) { fUseTailCleanup=val; }
00542   void SetNjet25(int val);
00543   void SetAllowUseHT(bool allowUseHT) { rawInput.allowUseHT=allowUseHT;}
00544   void SetMETScanScheme(int val) { METScanScheme=val; }
00545 
00546   void SetMHtSigma1(double val) { rawInput.MHtSigma1=val; }
00547   void SetMHtSigma2(double val) { rawInput.MHtSigma2=val; } 
00548   void SetMHtGaussFr(double val) { rawInput.MHtGaussFr=val; }
00549   void SetMnuScanRange(double val) { MnuScanRange=val; }
00550 
00551   void SetProposalTryMEt(double val) {m_ProposalTryMEt=val; }
00552   void SetProposalTryPhi(double val) {m_ProposalTryPhi=val;}
00553   void SetProposalTryMnu(double val) {m_ProposalTryMnu=val;}
00554   void SetProposalTryEtau(double val) {m_ProposalTryEtau=val;}
00555 
00556   // added by Tomas Davidek
00557   void SetUseDphiLL(bool val) {fUseDphiLL = val;}
00558 
00559   void SetLFVmode(int val) { LFVmode=val; }
00560   MMCCalibrationSet::e GetCalibrationSet() { return m_MMCCalibrationSet; }
00561   int GetNiterFit1() const { return Niter_fit1; } // number of iterations per loop in dPhi loop
00562   int GetNiterFit2() const { return Niter_fit2; } // number of iterations per loop in MET loop
00563   int GetNiterFit3() const { return Niter_fit3; } // number of iterations per loop in Mnu loop
00564   int GetNiterRandom() const { return m_NiterRandomLocal; } // number of random iterations
00565 
00566   int GetNsucStop() const { return NsucStop; } // Arrest criteria for NSuc
00567   int GetRMSStop() const { return RMSStop; }
00568   double GetMeanbinStop() const { return m_meanbinStop;}
00569   int GetRndmSeedAltering() const { return RndmSeedAltering; } // number of iterations per loop in Mnu loop
00570 
00571   bool GetUseHT() const { return OutputInfo.UseHT; } // if use HT
00572 
00573   int GetMarkovCountDuplicate() const { return m_markovCountDuplicate; }
00574   int GetMarkovNRejectNoSol() const { return m_markovNRejectNoSol;}
00575   int GetMarkovNRejectMetropolis() const {return m_markovNRejectMetropolis;}
00576   int GetMarkovNAccept() const { return m_markovNAccept; }
00577   int GetMarkovNFullscan() const { return m_markovNFullScan;}
00578   double GetProposalTryMEt() const {return m_ProposalTryMEt;}
00579   double GetProposalTryPhi() const {return m_ProposalTryPhi;}
00580   double GetProposalTryMnu() const {return m_ProposalTryMnu;}
00581   double GetProposalTryEtau() const {return m_ProposalTryEtau;}
00582 
00583 
00584   double GetmMaxError() const {return m_PrintmMaxError;}
00585   double GetmMeanError() const { return m_PrintmMeanError;}
00586   double GetmInvWidth2Error() const {return m_PrintmInvWidth2Error;}
00587 
00588 
00589   //-------- Get results;
00590   int GetFitStatus(); // return fit status
00591   double GetFittedMass(int fitcode); // returns fitted Mass
00592   std::shared_ptr<TH1F> GetMassHistogram(); // return mass histogram
00593   std::shared_ptr<TH1F> GetMassHistogramNoWeight(); // return mass histogram without weights
00594   int GetNTrials();// total number of point scanned
00595   int GetNSuccesses();// total number of point with at least 1 solutions
00596   int GetNSolutions();// total number of solutions 
00597   double GetSumW(); // sum of weights
00598   double GetAveSolRMS(); // ave RMS of solutions (for one event)
00599 
00600   double GetRms2Mpv(); // returns RMS/MPV according to histogram method
00601   TLorentzVector GetNeutrino4vec(int fitcode, int ind); // returns neutrino 4-vec
00602   double GetFitSignificance(int fitcode); // returns fit significance
00603   TLorentzVector GetTau4vec(int fitcode, int ind); // returns full tau 4-vec
00604   TLorentzVector GetResonanceVec(int fitcode); // returns 4-vec for resonance 
00605   TVector2 GetFittedMetVec(int fitcode); // returns 2-vec for fitted MET
00606   const AllPointsOutputInfo & GetAllPointsOutputInfo() { return m_aPOI;};
00607   int StandardCollApprox(const TLorentzVector & tau_vec1, const TLorentzVector & tau_vec2, const TVector2 & met_vec, double &Mrec); // standard collinear approximation
00608   // compute maximum from histo
00609   //  double maxFromHist(TH1F * theHist, double & prob, double & h_integral, double & h_chi2, double & h_discri, double & h_tantheta, double & h_tanthetaw, double & h_fitlength, double & h_rms, double & h_rmsvsdiscri, const MaxHistStrategy::e maxHistStrategy=MaxHistStrategy::FIT,const int winHalfWidth=2,bool debug=false);
00610   Double_t maxFitting(Double_t *x, Double_t *par);
00611   
00612   double maxFromHist(TH1F *theHist, std::vector<double> & histInfo, const MaxHistStrategy::e maxHistStrategy=MaxHistStrategy::FIT,const int winHalfWidth=2,bool debug=false);
00613   double maxFromHist(const std::shared_ptr<TH1F>& theHist, std::vector<double> & histInfo, const MaxHistStrategy::e maxHistStrategy=MaxHistStrategy::FIT,const int winHalfWidth=2,bool debug=false) {
00614     return maxFromHist(theHist.get(), histInfo, maxHistStrategy, winHalfWidth, debug);
00615   }
00616 
00617     // missing removed interfaces to turn on LFV code
00618     TVector2 metvec_tmp;
00619     inline double dTheta3DLimit(const int & tau_type, const int & limit_code,const double & P_tau);
00620 
00621 };
00622 
00623 #endif