/////////////////////////////////////////////////////////////////////
//
// Resonance mass sampling from a relativistic
// Breit-Wigner distribution with mass-dependent width Breit-Wigner
// (low-mass cutoff through Gamma0*(q/q0)**(2L+1) ).
//
//                                  Author:  Kagarlis 
//                                  Implemented: I. Froehlich
//                                  Written: 23.5.2007
//                                  Revised: 
// 
/////////////////////////////////////////////////////////////////////


using namespace std;
#include <sstream>
#include <iostream>
#include <iomanip>

#include "PBreitWigner.h"

ClassImp(PBreitWigner)

PBreitWigner::PBreitWigner()  {
};

PBreitWigner::PBreitWigner(const Char_t *id, const Char_t *de, Int_t key) :
    PHadronModel(id, de, key) {

    width_model = NULL;
    mr = makeStaticData()->GetParticleMassByKey(primary_key);
};

PDistribution *PBreitWigner::Clone(const char *) const {
    return new PBreitWigner((const PBreitWigner &)* this);
};

Bool_t PBreitWigner::SampleMass(Double_t *mass, Int_t *didx) {
    // Resonance mass sampling from a relativistic
    // Breit-Wigner distribution with mass-dependent width Breit-Wigner
    //(low-mass cutoff through Gamma0*(q/q0)**(2L+1) ).
    // Used e.g. for 3-body decays in PHadronDecayM3.

    if (didx) {
	if (didx_option != didx[0]) {
	    didx_option = didx[0];
	    SetParameter(1, didx_option);
	}
    } else {
	if (didx_option != -1) {
	    didx_option = -1;
	    SetParameter(1, didx_option);
	}
    }
    
    mass[0] = this->GetRandom();

    return kTRUE;
}

Double_t PBreitWigner::GetWeight(Double_t *mass, Int_t *didx) {
// relativistic Breit-Wigner distribution function for particle "key"
    
    double m = mass[0];

    if ((m < GetMin()) || (m > GetMax())) return 0.;

    double m2 = m*m, 	
	mm = mr*mr-m2;
    int didx_local = -1;
    if (didx) didx_local = didx[0];

    global_weight_scaling = makeDynamicData()->GetParticleScalingFactor(is_pid);

    double wmt = 1.;
    if (!GetWidth(m, &width)) {
	Warning("GetWeight", "GetWidth failed");
	return -1;
    }
    //N.B.: If didx is used, we sample the partial decay width

    double g = width, 
	g2 = g*g;
    double partial_width;
    if (didx_local >= 0) {
	if (!GetWidth(m, &partial_width, didx_local)) {
	    Warning("GetWeight", "GetWidth failed");
	    return -1;
	}	
	if (width_model) {
	    width_model->GetWidth(m, &partial_width, didx_local);
	}
    } else 
	partial_width = width;


// I copied this comment from PData here for bookkeeping (IF)
//  if (tempMtScaling>0. && (id==41 || id==42 || id==43)) { // obsolete > v4.09
//    fold rho line shape with thermal distribution
//    double mcut = TMath::Max(m,0.28);  // cut mass at 2 pion masses
//    wmt = mtIntegral(mcut,tempMtScaling)/mtIntegral(Mass(id),tempMtScaling);
//  } 
//  return wmt*m2*g2/(mm*mm+m2*g2);  // ->  BW(Mres) = 1
//    cout << "result: " << global_weight_scaling*wmt*m2*partial_width/(mm*mm+m2*g2) << endl;

    Double_t w = global_weight_scaling*wmt*m2*partial_width/(mm*mm+m2*g2);

    return w;
}

 PBreitWigner.cc:1
 PBreitWigner.cc:2
 PBreitWigner.cc:3
 PBreitWigner.cc:4
 PBreitWigner.cc:5
 PBreitWigner.cc:6
 PBreitWigner.cc:7
 PBreitWigner.cc:8
 PBreitWigner.cc:9
 PBreitWigner.cc:10
 PBreitWigner.cc:11
 PBreitWigner.cc:12
 PBreitWigner.cc:13
 PBreitWigner.cc:14
 PBreitWigner.cc:15
 PBreitWigner.cc:16
 PBreitWigner.cc:17
 PBreitWigner.cc:18
 PBreitWigner.cc:19
 PBreitWigner.cc:20
 PBreitWigner.cc:21
 PBreitWigner.cc:22
 PBreitWigner.cc:23
 PBreitWigner.cc:24
 PBreitWigner.cc:25
 PBreitWigner.cc:26
 PBreitWigner.cc:27
 PBreitWigner.cc:28
 PBreitWigner.cc:29
 PBreitWigner.cc:30
 PBreitWigner.cc:31
 PBreitWigner.cc:32
 PBreitWigner.cc:33
 PBreitWigner.cc:34
 PBreitWigner.cc:35
 PBreitWigner.cc:36
 PBreitWigner.cc:37
 PBreitWigner.cc:38
 PBreitWigner.cc:39
 PBreitWigner.cc:40
 PBreitWigner.cc:41
 PBreitWigner.cc:42
 PBreitWigner.cc:43
 PBreitWigner.cc:44
 PBreitWigner.cc:45
 PBreitWigner.cc:46
 PBreitWigner.cc:47
 PBreitWigner.cc:48
 PBreitWigner.cc:49
 PBreitWigner.cc:50
 PBreitWigner.cc:51
 PBreitWigner.cc:52
 PBreitWigner.cc:53
 PBreitWigner.cc:54
 PBreitWigner.cc:55
 PBreitWigner.cc:56
 PBreitWigner.cc:57
 PBreitWigner.cc:58
 PBreitWigner.cc:59
 PBreitWigner.cc:60
 PBreitWigner.cc:61
 PBreitWigner.cc:62
 PBreitWigner.cc:63
 PBreitWigner.cc:64
 PBreitWigner.cc:65
 PBreitWigner.cc:66
 PBreitWigner.cc:67
 PBreitWigner.cc:68
 PBreitWigner.cc:69
 PBreitWigner.cc:70
 PBreitWigner.cc:71
 PBreitWigner.cc:72
 PBreitWigner.cc:73
 PBreitWigner.cc:74
 PBreitWigner.cc:75
 PBreitWigner.cc:76
 PBreitWigner.cc:77
 PBreitWigner.cc:78
 PBreitWigner.cc:79
 PBreitWigner.cc:80
 PBreitWigner.cc:81
 PBreitWigner.cc:82
 PBreitWigner.cc:83
 PBreitWigner.cc:84
 PBreitWigner.cc:85
 PBreitWigner.cc:86
 PBreitWigner.cc:87
 PBreitWigner.cc:88
 PBreitWigner.cc:89
 PBreitWigner.cc:90
 PBreitWigner.cc:91
 PBreitWigner.cc:92
 PBreitWigner.cc:93
 PBreitWigner.cc:94
 PBreitWigner.cc:95
 PBreitWigner.cc:96
 PBreitWigner.cc:97
 PBreitWigner.cc:98
 PBreitWigner.cc:99
 PBreitWigner.cc:100
 PBreitWigner.cc:101
 PBreitWigner.cc:102
 PBreitWigner.cc:103
 PBreitWigner.cc:104
 PBreitWigner.cc:105
 PBreitWigner.cc:106
 PBreitWigner.cc:107
 PBreitWigner.cc:108
 PBreitWigner.cc:109
 PBreitWigner.cc:110