RunningWidthCalculator Class Reference

#include <RunningWidthCalculator.h>

Public Member Functions
	RunningWidthCalculator (const DalitzEventPattern &pat)
TH1D *	makeHisto_phaseSpace (int nBins, double max_m_inMeV, std::string OutputFileName="")
TH1D *	makeHisto_2body (int nBins, double max_m_inMeV, const DecayTree &thisDcy, std::string OutputFileName="")
TH1D *	makeHisto_3body (int nBins, double max_m_inMeV, const DecayTree &thisDcy, std::string OutputFileName="")
TH1D *	makeHisto_dalitz (int nBins, double max_m_inMeV, int nIntegrationEvents, IFastAmplitudeIntegrable *amps=0, std::string OutputFileName="")
TF1 *	getRunningWidthFunction_phaseSpace (double max_s_inGeV2)
TF1 *	getRunningWidthFunction_2body (double max_s_inGeV2, const DecayTree &thisDcy)
TF1 *	getRunningWidthFunction_3body (double max_s_inGeV2, const DecayTree &thisDcy)
TF1 *	getRunningWidthFunction_Fermi (double max_s_inGeV2, double lambda, double s0)
std::vector< double >	getPartialWidthCouplingsFromBF (std::vector< double > BF, std::vector< TF1 * > partialWidths, double m0, double gamma0)
TH1D *	makeRunningMassHisto_3body (int nBins, double max_s_inGeV2, std::vector< TF1 * > gamma, std::vector< double > g, std::string OutputFileName="")
void	setDalitzEventPattern (const DalitzEventPattern &pat)

Private Member Functions
void	set_min_s_inGeV2 ()

Private Attributes
DalitzEventPattern	_pat
double	_min_s_inGeV2

Detailed Description

Definition at line 13 of file RunningWidthCalculator.h.

Constructor & Destructor Documentation

RunningWidthCalculator()

RunningWidthCalculator::RunningWidthCalculator

(

const DalitzEventPattern &

pat

)

Definition at line 37 of file RunningWidthCalculator.cpp.

                                                                           :
    _pat(pat)
{
  cout << "RunningWidthCalculator::RunningWidthCalculator: I got called." << endl;
    set_min_s_inGeV2();
}

Member Function Documentation

getPartialWidthCouplingsFromBF()

std::vector< double > RunningWidthCalculator::getPartialWidthCouplingsFromBF	(	std::vector< double >	BF,
		std::vector< TF1 * >	partialWidths,
		double	m0,
		double	gamma0
	)

Definition at line 452 of file RunningWidthCalculator.cpp.

                                                                                                                                                      {

    //MinuitParameterSet mps();
    
    std::vector<FitParameter*> fit_couplings;
    
    FitParameter* scale = new FitParameter("scale",2,1,0.1);
    fit_couplings.push_back(scale);
    
    for (unsigned int i=1; i < BF.size() ; i++) {
        FitParameter* g =
    new FitParameter(("g_"+anythingToString((int) i)).c_str(),0,1,0.1);
    //new FitParameter(("g_"+std::to_string(i)).c_str(),0,1,0.1);
        fit_couplings.push_back(g);
    }
      
    chi2BF f(BF, partialWidths, fit_couplings, m0, gamma0);
    Minimiser mini(&f);
    mini.doFit();
    
    vector<double> couplings;
    double norm = 0.;
    
    
    double sum = 0.;
    for (unsigned int i=1; i<fit_couplings.size(); i++) {
        sum += *fit_couplings[i];
    }
    
    for (unsigned int i=0; i<fit_couplings.size(); i++) {
        double c;
        if(i==0) c = 1.;//(1.-sum);
        else c =  *fit_couplings[i];
        norm += c*partialWidths[i]->Eval(m0*m0/(GeV*GeV));
    }
    
    for (unsigned int i=0; i<fit_couplings.size(); i++) {
        double c;
        if(i==0) c = 1;//(1.-sum);
        else c =  *fit_couplings[i];
        couplings.push_back(c / norm * gamma0/GeV);
    }
    
    return couplings;
     
}

getRunningWidthFunction_2body()

TF1 * RunningWidthCalculator::getRunningWidthFunction_2body	(	double	max_s_inGeV2,
		const DecayTree &	thisDcy
	)

Definition at line 564 of file RunningWidthCalculator.cpp.

                                                                                                       {
    
    bool dbThis=false;
    
    // Some safety checks
    if(thisDcy.finalState().size() != 2){
        cout << "ERROR in RunningWidthCalculator::makeHisto_2body: I can handle only 2 body decays but the decay " << thisDcy 
        << " has " << thisDcy.finalState().size() << " final state particles. I'll return 0." << endl;
        return 0;
    }
    
    if(thisDcy.nDgtr() != 2){
        cout << "ERROR in RunningWidthCalculator::makeHisto_2body: "
        << " Mother should have 2 daughters, but it says it has "
        << thisDcy.nDgtr() << ". I'll return 0." << endl;
        return 0;
    }
    
    if( ! _pat.compatibleWith(thisDcy)){
        cout << "ERROR in RunningWidthCalculator::makeHisto_2body:"
        << "DalitzEventPattern and DecayTree are not compatible. I'll return 0." << endl;
        return 0;
    }
        
    // Use the BW_BW class to add resonances to ResonanceProperitiesList 
    // and to get acces to angular momentum 
    AssociatingDecayTree associatingDecayTreeX(thisDcy);
    const AssociatedDecayTree X = associatingDecayTreeX.getTree(_pat);
    BW_BW BW_X(X);
    
    // Define all parameters
    // 2 body decay : X -> 1 2
    const int nPar = 4;
    Double_t par[nPar];
    
    // Final state (1,2) masses:   
    par[0] = ParticlePropertiesList::getMe()->get((int) X.getDgtrTreePtr(0)->getVal())->mass();
    par[1] = ParticlePropertiesList::getMe()->get((int) X.getDgtrTreePtr(1)->getVal())->mass();
    // Angular momentum in X -> 1 2
    par[2] = (BW_X.twoLPlusOne()-1)/2; 
    // Mother (X) radius
    par[3] = ResonancePropertiesList::getMe()->get((int)_pat[0])->radius();
        
    TF1 *gamma = new TF1("Gamma_2body",Gamma_2body,0,max_s_inGeV2,nPar);
    gamma->SetParameters(par); 
    
    return gamma;
    (void)dbThis;
}

getRunningWidthFunction_3body()

TF1 * RunningWidthCalculator::getRunningWidthFunction_3body	(	double	max_s_inGeV2,
		const DecayTree &	thisDcy
	)

Definition at line 285 of file RunningWidthCalculator.cpp.

                                                                                                       {

    bool dbThis=false;
    
    // Some safety checks
    if(thisDcy.finalState().size() != 3){
        cout << "ERROR in RunningWidthCalculator::makeHisto_3body: I can handle only 3 body decays but the decay " << thisDcy 
        << " has " << thisDcy.finalState().size() << " final state particles. I'll return 0." << endl;
        return 0;
    }
    
    if(thisDcy.nDgtr() != 2){
        cout << "ERROR in RunningWidthCalculator::makeHisto_3body: "
        << " Mother should have 2 daughters, but it says it has "
        << thisDcy.nDgtr() << ". I'll return 0." << endl;
        return 0;
    }
    
    if( ! _pat.compatibleWith(thisDcy)){
        cout << "ERROR in RunningWidthCalculator::makeHisto_3body:"
        << "DalitzEventPattern and DecayTree are not compatible. I'll return 0." << endl;
        return 0;
    }
    
    // Everything is fine, now search the resonance
    MINT::const_counted_ptr<AssociatedDecayTree> R, bachelor;
    for(int i=0; i< thisDcy.nDgtr(); i++){
        const_counted_ptr<AssociatedDecayTree> dgtr = thisDcy.getDgtrTreePtr(i);
        if(dgtr->nDgtr() == 2) R = dgtr;
        else bachelor = dgtr; 
    }
    
    if(0==R || 0== bachelor){
        cout << "ERROR in RunningWidthCalculator::makeHisto_3body:"
        << " Didn't find resonance or bachelor. I'll return 0.  " << endl;
        return 0;
    }
    
    // Use the BW_BW class to add resonances to ResonanceProperitiesList 
    // and to get acces to angular momentum 
    AssociatingDecayTree associatingDecayTreeX(thisDcy);
    const AssociatedDecayTree X = associatingDecayTreeX.getTree(_pat);
    
    BW_BW BW_X(X);
    BW_BW BW_R(*R);
    
    // Define all parameters
    // 3 body decay : X -> (R->1 2) 3
    const int nPar = 10;
    Double_t par[nPar];
    
    // Mother (X) radius
    par[0] = ResonancePropertiesList::getMe()->get((int)_pat[0])->radius();
    // Final state (1,2,3) masses:   
    par[1] = ParticlePropertiesList::getMe()->get((int) R->getDgtrTreePtr(0)->getVal())->mass();
    par[2] = ParticlePropertiesList::getMe()->get((int) R->getDgtrTreePtr(1)->getVal())->mass();
    par[3] = ParticlePropertiesList::getMe()->get((int) bachelor->getVal())->mass();
    // Angular momentum in X -> R 3
    par[4] = (BW_X.twoLPlusOne()-1)/2; 
    // Angular momentum in R -> 1 2
    par[5] = (BW_R.twoLPlusOne()-1)/2; 
    // Resonance (R) mass, width and radius
    par[6] = ResonancePropertiesList::getMe()->get((int) R->getVal())->mass();
    par[7] = ResonancePropertiesList::getMe()->get((int) R->getVal())->width();
    par[8] = ResonancePropertiesList::getMe()->get((int) R->getVal())->radius();
    // Upper limit of s(1,2,3)
    par[9] = max_s_inGeV2 * GeV *GeV;
    
    if(dbThis) for(int i=0; i<nPar; i++) cout << "Parameter " << i << " = " << par[i] << endl;
    
    TF1 *gamma = new TF1("Gamma_mother_3body",Gamma_mother_3body,_min_s_inGeV2,max_s_inGeV2,nPar);
    gamma->SetParameters(par); 
    
    return gamma;
}

getRunningWidthFunction_Fermi()

TF1 * RunningWidthCalculator::getRunningWidthFunction_Fermi	(	double	max_s_inGeV2,
		double	lambda,
		double	s0
	)

Definition at line 758 of file RunningWidthCalculator.cpp.

                                                                                                       {
    
    // Define all parameters
    const int nPar = 3;
    Double_t par[nPar];
    par[0] = lambda;
    par[1] = s0;
    par[2] = _min_s_inGeV2;

    TF1 *gamma = new TF1("Fermi_phaseSpace",Fermi_phaseSpace,0,max_s_inGeV2,nPar);
    gamma->SetParameters(par); 
    
    return gamma;
}

getRunningWidthFunction_phaseSpace()

TF1 * RunningWidthCalculator::getRunningWidthFunction_phaseSpace

(

double

max_s_inGeV2

)

Definition at line 728 of file RunningWidthCalculator.cpp.

                                                                                  {
    
    if(_pat.size() != 4) {
        cout << "ERROR in RunningWidthCalculator::makeHisto_phaseSpace: I can handle only 3 body decays but the pattern " << _pat 
        << " has " << _pat.size()-1 << " final state particles. I'll return 0." << endl;
        return 0;    
    }
    
    const int nPar = 4;
    Double_t par[nPar];
    par[0] = ParticlePropertiesList::getMe()->get((int)_pat[0])->mass();
    par[1] = ParticlePropertiesList::getMe()->get((int)_pat[1])->mass();
    par[2] = ParticlePropertiesList::getMe()->get((int)_pat[2])->mass();
    par[3] = ParticlePropertiesList::getMe()->get((int)_pat[3])->mass();
    
    TF1 *ps = new TF1("phaseSpaceIntegral",phaseSpaceIntegral,0,max_s_inGeV2,nPar);
    ps->SetParameters(par);
    
    return ps;    
}

makeHisto_2body()

TH1D * RunningWidthCalculator::makeHisto_2body	(	int	nBins,
		double	max_m_inMeV,
		const DecayTree &	thisDcy,
		std::string	OutputFileName = ""
	)

Definition at line 615 of file RunningWidthCalculator.cpp.

                                                                                                                            {
    
    bool dbThis=false;
    
    // Some safety checks
    if(thisDcy.finalState().size() != 2){
        cout << "ERROR in RunningWidthCalculator::makeHisto_2body: I can handle only 2 body decays but the decay " << thisDcy 
        << " has " << thisDcy.finalState().size() << " final state particles. I'll return 0." << endl;
        return 0;
    }
    
    if(thisDcy.nDgtr() != 2){
        cout << "ERROR in RunningWidthCalculator::makeHisto_2body: "
        << " Mother should have 2 daughters, but it says it has "
        << thisDcy.nDgtr() << ". I'll return 0." << endl;
        return 0;
    }
    
    if( ! _pat.compatibleWith(thisDcy)){
        cout << "ERROR in RunningWidthCalculator::makeHisto_2body:"
        << "DalitzEventPattern and DecayTree are not compatible. I'll return 0." << endl;
        return 0;
    }
    
    // Use the BW_BW class to add resonances to ResonanceProperitiesList 
    // and to get acces to angular momentum 
    AssociatingDecayTree associatingDecayTreeX(thisDcy);
    const AssociatedDecayTree X = associatingDecayTreeX.getTree(_pat);
    BW_BW BW_X(X);
    
    // Define all parameters
    // 2 body decay : X -> 1 2
    //const int nPar = 10;
    //Double_t par[nPar];
    
    // Mother (X) radius
    double r = ResonancePropertiesList::getMe()->get((int)_pat[0])->radius();
    // Final state (1,2) masses:   
    double m1 = ParticlePropertiesList::getMe()->get((int) X.getDgtrTreePtr(0)->getVal())->mass();
    double m2 = ParticlePropertiesList::getMe()->get((int) X.getDgtrTreePtr(1)->getVal())->mass();
    // Angular momentum in X -> 1 2
    double l = (BW_X.twoLPlusOne()-1)/2; 
        
    if(OutputFileName=="")
        OutputFileName = "RunningWidth_"+ParticlePropertiesList::getMe()->get(abs((int)_pat[0]))->name()+"_2body.root";
    TFile* f= new TFile(OutputFileName.c_str(),"RECREATE");  
    
    TH1D* h= new TH1D("RunningWidth",";s [GeV^2]; #Gamma (s) [a.u.]",nBins,_min_s_inGeV2,max_s_inGeV2); 
    for (int i=1; i<= h->GetNbinsX(); i++) {
        h->SetBinContent(i, Gamma_2body( sqrt(h->GetXaxis()->GetBinCenter(i))*GeV , m1,m2, l, r));
    }  
    
    //h->Scale(1./h->Integral());    
    h->Write();
    f->Write();  
    return h;
    (void)dbThis;
}

makeHisto_3body()

TH1D * RunningWidthCalculator::makeHisto_3body	(	int	nBins,
		double	max_m_inMeV,
		const DecayTree &	thisDcy,
		std::string	OutputFileName = ""
	)

Definition at line 196 of file RunningWidthCalculator.cpp.

                                                                                                                            {
     
     bool dbThis=false;
    
     // Some safety checks
     if(thisDcy.finalState().size() != 3){
         cout << "ERROR in RunningWidthCalculator::makeHisto_3body: I can handle only 3 body decays but the decay " << thisDcy 
         << " has " << thisDcy.finalState().size() << " final state particles. I'll return 0." << endl;
         return 0;
     }
    
     if(thisDcy.nDgtr() != 2){
        cout << "ERROR in RunningWidthCalculator::makeHisto_3body: "
        << " Mother should have 2 daughters, but it says it has "
        << thisDcy.nDgtr() << ". I'll return 0." << endl;
        return 0;
     }
     
     if( ! _pat.compatibleWith(thisDcy)){
         cout << "ERROR in RunningWidthCalculator::makeHisto_3body:"
         << "DalitzEventPattern and DecayTree are not compatible. I'll return 0." << endl;
         return 0;
     }
    
     // Everything is fine, now search the resonance
     MINT::const_counted_ptr<AssociatedDecayTree> R, bachelor;
     for(int i=0; i< thisDcy.nDgtr(); i++){
        const_counted_ptr<AssociatedDecayTree> dgtr = thisDcy.getDgtrTreePtr(i);
        if(dgtr->nDgtr() == 2) R = dgtr;
        else bachelor = dgtr; 
     }
    
     if(0==R || 0== bachelor){
        cout << "ERROR in RunningWidthCalculator::makeHisto_3body:"
        << " Didn't find resonance or bachelor. I'll return 0.  " << endl;
        return 0;
     }
    
     // Use the BW_BW class to add resonances to ResonanceProperitiesList 
     // and to get acces to angular momentum 
     AssociatingDecayTree associatingDecayTreeX(thisDcy);
     const AssociatedDecayTree X = associatingDecayTreeX.getTree(_pat);
    
     BW_BW BW_X(X);
     BW_BW BW_R(*R);
    
     // Define all parameters
     // 3 body decay : X -> (R->1 2) 3
     const int nPar = 10;
     Double_t par[nPar];
    
     // Mother (X) radius
     par[0] = ResonancePropertiesList::getMe()->get((int)_pat[0])->radius();
     // Final state (1,2,3) masses:   
     par[1] = ParticlePropertiesList::getMe()->get((int) R->getDgtrTreePtr(0)->getVal())->mass();
     par[2] = ParticlePropertiesList::getMe()->get((int) R->getDgtrTreePtr(1)->getVal())->mass();
     par[3] = ParticlePropertiesList::getMe()->get((int) bachelor->getVal())->mass();
     // Angular momentum in X -> R 3
     par[4] = (BW_X.twoLPlusOne()-1)/2; 
     // Angular momentum in R -> 1 2
     par[5] = (BW_R.twoLPlusOne()-1)/2; 
     // Resonance (R) mass, width and radius
     par[6] = ResonancePropertiesList::getMe()->get((int) R->getVal())->mass();
     par[7] = ResonancePropertiesList::getMe()->get((int) R->getVal())->width();
     par[8] = ResonancePropertiesList::getMe()->get((int) R->getVal())->radius();
     // Upper limit of s(1,2,3)
     par[9] = max_s_inGeV2 * GeV *GeV;
    
     if(dbThis) for(int i=0; i<nPar; i++) cout << "Parameter " << i << " = " << par[i] << endl;
     
    if(OutputFileName=="")
        OutputFileName = "RunningWidth_"+ParticlePropertiesList::getMe()->get(abs((int)_pat[0]))->name()+"_3body.root";
     TFile* f= new TFile(OutputFileName.c_str(),"RECREATE");  
     
     TF1 *gamma = new TF1("Gamma_mother_3body",Gamma_mother_3body,0,max_s_inGeV2,nPar);
     gamma->SetParameters(par); 
     
     TH1D* h= new TH1D("RunningWidth",";s [GeV^2]; #Gamma (s) [a.u.]",nBins,0,max_s_inGeV2); 
     for (int i=1; i<= h->GetNbinsX(); i++) {
         h->SetBinContent(i,gamma->Eval(h->GetXaxis()->GetBinCenter(i) ));    
     }  
     
     //h->Scale(1./h->Integral());    
     h->Write();
     f->Write();  
     return h;
}

makeHisto_dalitz()

TH1D * RunningWidthCalculator::makeHisto_dalitz	(	int	nBins,
		double	max_m_inMeV,
		int	nIntegrationEvents,
		IFastAmplitudeIntegrable *	amps = 0,
		std::string	OutputFileName = ""
	)

Definition at line 50 of file RunningWidthCalculator.cpp.

                                                                                                                                                           {
    
    if(_pat.size() != 4) {
        cout << "ERROR in RunningWidthCalculator::makeHisto_dalitz: I can handle only 3 body decays but the pattern " << _pat 
        << " has " << _pat.size()-1 << " final state particles. I'll return 0." << endl;
        return 0;    
    }
    
    if(amps==0) amps = new FitAmpSum(_pat);
    //Important! Ensures everything is initialised
    DalitzEventList eventTest;
    eventTest.generatePhaseSpaceEvents(1,_pat);
    amps->RealVal(eventTest[0]);  
    
    if(OutputFileName=="")
        OutputFileName = "RunningWidth_"+ParticlePropertiesList::getMe()->get(abs((int)_pat[0]))->name()+"_Dalitz.root";
    
    TFile* f= new TFile(OutputFileName.c_str(),"RECREATE");     
    TH1D* h= new TH1D("RunningWidth","RunningWidth",nBins,_min_s_inGeV2,max_s_inGeV2); 
      
    const ParticleProperties* pp = ParticlePropertiesList::getMe()->get((int)_pat[0]);
    const ResonanceProperties* rp = ResonancePropertiesList::getMe()->get((int)_pat[0]);
        
    for (int b=1; b<= h->GetNbinsX(); b++) {  
        // Set mass to new value in ResonancePropertiesList
        // and in ParticlePropertiesList !
        // This is necessary since SignalGenerator uses the mass in ParticlePropertiesList
        // but it should use the mass in ResonancePropertiesList. Fix some day.
        rp->changeMassForDebug(sqrt(h->GetXaxis()->GetBinCenter(b))*GeV);
        pp->setMass(sqrt(h->GetXaxis()->GetBinCenter(b))*GeV);
        
        // Generate integration events
        SignalGenerator sg(_pat,amps);
        sg.setWeighted();
        sg.noPrintout();
        DalitzEventList eventList;
        // Maybe we should somehow scale the number of integration events with mass, 
        // i.e. the available phase space ?
        sg.FillEventList(eventList, nIntegrationEvents);

        // Calculate the integral
        double integral = 0.;
        double sumw = 0.;
        for (unsigned int i= 0; i< eventList.size(); i++) {
            integral += amps->RealVal(eventList[i])*eventList[i].getWeight()/(eventList[i].getGeneratorPdfRelativeToPhaseSpace());
            sumw += eventList[i].getWeight()/(eventList[i].getGeneratorPdfRelativeToPhaseSpace());
        }
        
        PhaseSpaceIntegral3body ps;
        double phaseSpace = ps.getVal(_pat);
                
        h->SetBinContent(b,integral*phaseSpace/sumw);    
    }

    h->Scale(h->GetNbinsX()/h->Integral());
    
    h->Write();
    f->Write();  
    
    return h;
}

makeHisto_phaseSpace()

TH1D * RunningWidthCalculator::makeHisto_phaseSpace	(	int	nBins,
		double	max_m_inMeV,
		std::string	OutputFileName = ""
	)

Definition at line 692 of file RunningWidthCalculator.cpp.

                                                                                                       {
    
    if(_pat.size() != 4) {
        cout << "ERROR in RunningWidthCalculator::makeHisto_phaseSpace: I can handle only 3 body decays but the pattern " << _pat 
        << " has " << _pat.size()-1 << " final state particles. I'll return 0." << endl;
        return 0;    
    }
    
    if(OutputFileName=="")
        OutputFileName = "RunningWidth_"+ParticlePropertiesList::getMe()->get(abs((int)_pat[0]))->name()+"_PhaseSpace.root";
    
    
    TFile* f= new TFile(OutputFileName.c_str(),"RECREATE");  
    TH1D* h= new TH1D("RunningWidth","RunningWidth",nBins,_min_s_inGeV2,max_s_inGeV2); 
        
    const int nPar = 4;
    Double_t par[nPar];
    par[0] = ParticlePropertiesList::getMe()->get((int)_pat[0])->mass();
    par[1] = ParticlePropertiesList::getMe()->get((int)_pat[1])->mass();
    par[2] = ParticlePropertiesList::getMe()->get((int)_pat[2])->mass();
    par[3] = ParticlePropertiesList::getMe()->get((int)_pat[3])->mass();
    
    TF1 *ps = new TF1("phaseSpaceIntegral",phaseSpaceIntegral,_min_s_inGeV2,max_s_inGeV2,nPar);
    ps->SetParameters(par);
    
    for (int i=1; i<= h->GetNbinsX(); i++) {
        h->SetBinContent(i,ps->Eval(h->GetXaxis()->GetBinCenter(i)));    
    }  
    
    h->Scale(h->GetNbinsX()/h->Integral());    
    h->Write();
    f->Write();  
    
    return h;    
}

makeRunningMassHisto_3body()

TH1D * RunningWidthCalculator::makeRunningMassHisto_3body	(	int	nBins,
		double	max_s_inGeV2,
		std::vector< TF1 * >	gamma,
		std::vector< double >	g,
		std::string	OutputFileName = ""
	)

Definition at line 515 of file RunningWidthCalculator.cpp.

                                                                                                                                                   {

    double absTol = 0.0001; 
    double relTol = 0.001; 
    unsigned int size = 1000; 
    int rule = 3;  
    double epsilon= 0.00001;
    
    if(OutputFileName=="")
        OutputFileName = "RunningMass_"+ParticlePropertiesList::getMe()->get(abs((int)_pat[0]))->name()+"_3body.root";
    TFile* f= new TFile(OutputFileName.c_str(),"RECREATE");  
    TH1D* h_m= new TH1D("RunningMass",";s [GeV^{2}]; m^{2} (s) [GeV^{2}]",nBins,_min_s_inGeV2,max_s_inGeV2); 
    for (int i=1; i<= h_m->GetNbinsX(); i++) {
        double singularity = h_m->GetXaxis()->GetBinCenter(i);
        if(singularity<_min_s_inGeV2){
            h_m->SetBinContent(i,0);
            continue;
        }
        Dispersion_Integrand integrand(gamma,g,singularity);
        ROOT::Math::IntegratorOneDim integrator(integrand, ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR,  absTol , relTol , size ,  rule ) ;
        h_m->SetBinContent(i,(integrator.Integral(_min_s_inGeV2,singularity-epsilon) + integrator.Integral(singularity+epsilon,max_s_inGeV2*10000))/TMath::Pi()); 
    }
    
    //h_m->Scale(h_m->GetNbinsX()/abs(h_m->Integral()));    
    h_m->Write();
    f->Write();  
    
    return h_m;
    /*
     double singularity = 1.4;
     
     Dispersion_Integrand integrand(gamma,singularity,100);
     ROOT::Math::IntegratorOneDim integrator(integrand, ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR,  absTol , relTol , size ,  rule )  ;
     double singularPts[3] = {_min_m_squared,singularity,singularity};
     std::vector<double> sp(singularPts, singularPts+3);  
     
     cout << "Int = " << integrator.Integral(_min_m_squared,singularity-epsilon) + integrator.IntegralUp(singularity+epsilon) << endl;
     
     
     Dispersion_Integrand_Subtraction integrand_sub(gamma,singularity);
     ROOT::Math::IntegratorOneDim integrator_sub(integrand_sub, ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR,  absTol , relTol , size ,  rule )  ;
     double tmp = integrator_sub.Integral(_min_m_squared,3);
     
     cout << "Int = " << -1.*(tmp + gamma->Eval(singularity)* log( (3-singularity)/(singularity-_min_m_squared) ) ) + integrator.IntegralUp(3)  << endl;
     */
}

set_min_s_inGeV2()

void RunningWidthCalculator::set_min_s_inGeV2 ( )

private

Definition at line 28 of file RunningWidthCalculator.cpp.

                                             {
    // Calculate lower phase space limit
    double min_m = 0.;
    for(unsigned int i= 1; i<_pat.size(); i++){
        min_m += _pat[i].mass()/GeV;
    }
    _min_s_inGeV2 = min_m*min_m;
}

setDalitzEventPattern()

void RunningWidthCalculator::setDalitzEventPattern

(

const DalitzEventPattern &

pat

)

Definition at line 44 of file RunningWidthCalculator.cpp.

                                                                               {
    _pat = pat;
    set_min_s_inGeV2();
}

Member Data Documentation

_min_s_inGeV2

double RunningWidthCalculator::_min_s_inGeV2

private

Definition at line 16 of file RunningWidthCalculator.h.

_pat

DalitzEventPattern RunningWidthCalculator::_pat

private

Definition at line 15 of file RunningWidthCalculator.h.

---

The documentation for this class was generated from the following files:

• Mint/RunningWidthCalculator.h
• src/Mojito/Lineshapes/RunningWidthCalculator.cpp