delete fMCtemplate;
}
//_________________________________________________________________________________________________
-Int_t AliAnalysisBtoJPSItoEle::DoMinimization(Double_t* x,
- Double_t* m, Int_t ncand)
+Int_t AliAnalysisBtoJPSItoEle::DoMinimization()
{
//
// performs the minimization
//
- AliInfo(Form("Number of candidates used for the minimisation is %d",ncand));
- fFCNfunction = new AliBtoJPSItoEleCDFfitHandler(x,m,ncand);
- SetResolutionConstants(fPtBin);
- SetCsiMC(fMCtemplate);
- fFCNfunction->SetErrorDef(0.5); // tells Minuit that the error interval is the one in which
- // the function differs from the minimum for less than setted value
Int_t iret=fFCNfunction->DoMinimization();
return iret;
}
//_________________________________________________________________________________________________
-void AliAnalysisBtoJPSItoEle::SetResolutionConstants(Int_t BinNum)
-{
- //
- // sets constants for parametrized resolution function
- //
- if(!fFCNfunction) {
- AliInfo("fFCNfunction not istanziated ---> nothing done");
- return;
- }
- AliInfo("Call likelihood SetResolutionConstants method ---> OK");
- AliBtoJPSItoEleCDFfitFCN* loglikePnt = fFCNfunction->LikelihoodPointer();
- if(!loglikePnt) {
- AliWarning("Pointer to AliBtoJPSItoEleCDFfitFCN class not found!");
- return;
- }
- loglikePnt->SetResolutionConstants(BinNum);
-}
-//_________________________________________________________________________________________________
void AliAnalysisBtoJPSItoEle::ReadCandidates(TNtuple* nt, Double_t* &pseudoproper,Double_t* &invmass, Int_t& ncand)
{
//
return;
}
//_________________________________________________________________________________________________
-void AliAnalysisBtoJPSItoEle::SetCsiMC(TH1F* MCtemplates)
+void AliAnalysisBtoJPSItoEle::SetCsiMC()
{
//
// Sets X distribution used as MC template for JPSI from B
//
- fFCNfunction->LikelihoodPointer()->SetCsiMC(MCtemplates);
+ fFCNfunction->LikelihoodPointer()->SetCsiMC(fMCtemplate);
return;
}
//_________________________________________________________________________________________________
+void AliAnalysisBtoJPSItoEle::SetFitHandler(Double_t* x /*pseudoproper*/, Double_t* m /*inv mass*/, Int_t ncand /*candidates*/)
+{
+ //
+ // Create the fit handler object to play with different params of the fitting function
+ //
+
+ fFCNfunction = new AliBtoJPSItoEleCDFfitHandler(x,m,ncand);
+ if(!fFCNfunction) {
+
+ AliInfo("fFCNfunction not istanziated ---> nothing done");
+ return;
+
+ }
+}
// Contact: Carmelo.Digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it\r
//-------------------------------------------------------------------------\r
class TNtuple ;\r
-class AliBtoJPSItoEleCDFfitFCN ;\r
-class AliBtoJPSItoEleCDFfitHandler ; \r
+//class AliBtoJPSItoEleCDFfitFCN ;\r
+//class AliBtoJPSItoEleCDFfitHandler ; \r
#include "TH1F.h"\r
+#include "AliBtoJPSItoEleCDFfitHandler.h"\r
+#include "AliBtoJPSItoEleCDFfitFCN.h"\r
\r
//_________________________________________________________________________ \r
class AliAnalysisBtoJPSItoEle : public TNamed {\r
AliAnalysisBtoJPSItoEle& operator=(const AliAnalysisBtoJPSItoEle& source);\r
virtual ~AliAnalysisBtoJPSItoEle();\r
\r
- Int_t DoMinimization(Double_t* x,Double_t* m, Int_t n);\r
+ Int_t DoMinimization();\r
void ReadCandidates(TNtuple* nt, Double_t* &x, Double_t* &m, Int_t& n); // primary JPSI + secondary JPSI + bkg couples\r
\r
void SetPtBin(Int_t BinNum) { fPtBin = BinNum ; }\r
- void SetCsiMC(TH1F* MCtemplate);\r
- void SetResolutionConstants(Int_t BinNum);\r
+ void SetCsiMC();\r
+ //void SetResolutionConstants(Int_t BinNum);\r
+ void SetFitHandler(Double_t* /*pseudoproper*/, Double_t* /*inv mass*/, Int_t /*candidates*/); \r
void CloneMCtemplate(const TH1F* MCtemplate) {fMCtemplate = (TH1F*)MCtemplate->Clone("fMCtemplate");}\r
Double_t* GetResolutionConstants(Double_t* resolutionConst);\r
+ AliBtoJPSItoEleCDFfitHandler* GetCDFFitHandler() const { return fFCNfunction ; }\r
Int_t GetPtBin() const { return fPtBin ; }\r
\r
private:\r
Int_t fPtBin; // number of pt bin in which the analysis is performes\r
TH1F* fMCtemplate; //! template of the MC distribution for the x distribution of the secondary J/psi\r
\r
- ClassDef(AliAnalysisBtoJPSItoEle,0); // AliAnalysisBtoJPSItoEle class\r
+ ClassDef(AliAnalysisBtoJPSItoEle,1); // AliAnalysisBtoJPSItoEle class\r
};\r
\r
#endif\r
void AliAnalysisTaskSEBtoJPSItoEle::Init()
{
// Initialization
- //Double_t ptCuts[2] = {1.,1.5}; // 1 GeV < pt < 1.5 GeV
- Double_t ptCuts[2] = {1.,100}; //
- SetPtCuts(ptCuts);
- SetMinimize(kTRUE);
- SetAODMCInfo(kTRUE);
if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::Init() \n");
fOutput->Add(fhDecayTimeMCjpsifromB);
}
- fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.);
+ fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",200,-2000.,2000.);
fhDecayTime->Sumw2();
fhDecayTime->SetMinimum(0);
fOutput->Add(fhDecayTime);
return;
}
+ // load MC particles
+ TClonesArray* mcArray =
+ dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
+ if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
+ AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));
+
// Read AOD Monte Carlo information
if(fOkAODMC) ReadAODMCInfo(aod,inputArrayJPSItoEle);
// loop over J/Psi->ee candidates
Int_t nInJPSItoEle = inputArrayJPSItoEle->GetEntriesFast();
- if(fDebug>1) printf("Number of B->JPSI->e+e-: %d\n",nInJPSItoEle);
+ printf("Number of JPSI->e+e- candidates: %d\n",nInJPSItoEle);
for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) {
AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayJPSItoEle->UncheckedAt(iJPSItoEle);
+
+ Int_t mcLabel = d->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
+
//Secondary vertex
Double_t vtxSec[3] = {0.,0.,0.};
Double_t vtxPrim[3] = {0.,0.,0.};
}
//Int_t okJPSI=0;
// here analyze only if cuts are passed
- //if(d->SelectBtoJPSI(fCuts,okJPSI)) {
+
if(d->Pt() < fPtCuts[1] && d->Pt() > fPtCuts[0]){ // apply pt cut only
+ if (mcLabel == -1)
+ {
+ AliDebug(2,"No MC particle found");
+ }
+ else {
+
fhInvMass->Fill(d->InvMassJPSIee());
fhD0->Fill(10000*d->ImpParXY());
fhD0D0->Fill(1e8*d->Prodd0d0());
fNtupleJPSI->Fill(pseudoProperDecayTime,d->InvMassJPSIee()); // fill N-tuple with X,M values
+ }//
+
} // end of JPSItoEle candidates selection according to cuts
if(unsetvtx) d->UnsetOwnPrimaryVtx();
}
printf("+++\n+++ MC template histo copied ---> OK\n+++\n");
- fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand);
+ //fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand);
+ fCdfFit->DoMinimization();
}
return;
for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i];
}
//_________________________________________________________________________________________________
-void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(const AliAODEvent* aodEvent, const TClonesArray* inputArray)
+void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray)
{
//
// Reads MC information if needed (for example in case of parametrized PID)
//
// load MC particles
- TClonesArray *mcArray =
- (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
- if(!mcArray) {
- printf("AliAnalysisTaskSEBtoJPSItoEle::UserExec: MC particles branch not found!\n");
- return;
- }
+ TClonesArray* mcArray =
+ dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+ if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
+ AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));
// load MC header
AliAODMCHeader* mcHeader =
printf("AliAnalysisTaskSEBtoJPSItoEle::UserExec: MC AOD header branch not found!\n");
}
- Double_t rmax = 0.000005;
+ AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
+
+ //Double_t rmax = 0.000005;
Double_t mcPrimVtx[3];
// get MC primary Vtx
Int_t labJPSIdaugh1=0;
for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) {
+
AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate);
+ Int_t mcLabelSec = dd->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
+
+ Double_t vtxPrim[3] = {0.,0.,0.};
+ Double_t vtxSec[3] = {0.,0.,0.};
+ dd->GetSecondaryVtx(vtxSec);
+ Bool_t unsetvtx=kFALSE;
+ if(!dd->GetOwnPrimaryVtx()) {
+ dd->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
+ unsetvtx=kTRUE;
+ }
+
if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only
AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0);
}
}
+ if (mcLabelSec == -1)
+ {
+ AliDebug(2,"No MC particle found");
+ }
+ else {
+
if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal
- AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0);
- pdgJPSI = partJPSI->GetPdgCode();
- if(pdgJPSI == 443){//this is for MC JPSI
- if(TMath::Sqrt((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Xv()-mcPrimVtx[0])+
- (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Yv()-mcPrimVtx[1])+
- (partJPSI->Zv()-mcPrimVtx[2])*(partJPSI->Zv()-mcPrimVtx[2])>rmax)){ //this is for MC displaced JPSI
- Int_t pdaugh0 = partJPSI->GetDaughter(0);
- Int_t pdaugh1 = partJPSI->GetDaughter(1);
- AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0);
- AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1);
- pdg0 = partDaugh0->GetPdgCode();
- pdg1 = partDaugh1->GetPdgCode();
- if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee
- labJPSIMother = partJPSI->GetMother();
- AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother);
- pdgJPSIMother = partJPSIMother->GetPdgCode();
- if(pdgJPSIMother==511 || pdgJPSIMother==521 ||
- pdgJPSIMother==10511 || pdgJPSIMother==10521 ||
- pdgJPSIMother==513 || pdgJPSIMother==523 ||
- pdgJPSIMother==10513 || pdgJPSIMother==10523 ||
- pdgJPSIMother==20513 || pdgJPSIMother==20523 ||
- pdgJPSIMother==515 || pdgJPSIMother==525 ||
- pdgJPSIMother==531 || pdgJPSIMother==10531 ||
- pdgJPSIMother==533 || pdgJPSIMother==10533 ||
- pdgJPSIMother==20533 || pdgJPSIMother==535 ||
- pdgJPSIMother==541 || pdgJPSIMother==10541 ||
- pdgJPSIMother==543 || pdgJPSIMother==10543 ||
- pdgJPSIMother==20543 || pdgJPSIMother==545)
- { //this is for MC JPSI -> ee from B-hadrons
+ AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0);
+ pdgJPSI = partJPSI->GetPdgCode();
+ Int_t pdaugh0 = partJPSI->GetDaughter(0);
+ Int_t pdaugh1 = partJPSI->GetDaughter(1);
+ AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0);
+ AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1);
+ pdg0 = partDaugh0->GetPdgCode();
+ pdg1 = partDaugh1->GetPdgCode();
+ if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee
+ labJPSIMother = partJPSI->GetMother();
+ AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother);
+ pdgJPSIMother = partJPSIMother->GetPdgCode();
+
+ // chech if JPSI comes from B
+ if( pdgJPSIMother==511 || pdgJPSIMother==521 ||
+ pdgJPSIMother==10511 || pdgJPSIMother==10521 ||
+ pdgJPSIMother==513 || pdgJPSIMother==523 ||
+ pdgJPSIMother==10513 || pdgJPSIMother==10523 ||
+ pdgJPSIMother==20513 || pdgJPSIMother==20523 ||
+ pdgJPSIMother==515 || pdgJPSIMother==525 ||
+ pdgJPSIMother==531 || pdgJPSIMother==10531 ||
+ pdgJPSIMother==533 || pdgJPSIMother==10533 ||
+ pdgJPSIMother==20533 || pdgJPSIMother==535 ||
+ pdgJPSIMother==541 || pdgJPSIMother==10541 ||
+ pdgJPSIMother==543 || pdgJPSIMother==10543 ||
+ pdgJPSIMother==20543 || pdgJPSIMother==545)
+ { //this is for MC JPSI -> ee from B-hadrons
+
+ fhInvMass->Fill(dd->InvMassJPSIee());
+ fhD0->Fill(10000*dd->ImpParXY());
+ fhD0D0->Fill(1e8*dd->Prodd0d0());
+ fhCosThetaStar->Fill(dd->CosThetaStarJPSI());
+ fhCtsVsD0D0->Fill(dd->CosThetaStarJPSI(),1e8*dd->Prodd0d0());
+ fhCosThetaPointing->Fill(dd->CosPointingAngle());
+
+ // compute X variable
+ AliAODVertex* primVertex = dd->GetOwnPrimaryVtx();
+ vtxPrim[0] = primVertex->GetX();
+ vtxPrim[1] = primVertex->GetY();
+ vtxPrim[2] = primVertex->GetZ();
+ Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dd->Px()) + (vtxSec[1]-vtxPrim[1])*(dd->Py()))/dd->Pt();
+ Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dd->Pt();
+
+ fhDecayTime->Fill(10000*pseudoProperDecayTime);
+ // Post the data already here
+ PostData(1,fOutput);
+
+ fNtupleJPSI->Fill(pseudoProperDecayTime,dd->InvMassJPSIee()); // fill N-tuple with X,M values
Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt();
Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt();
// Post the data already here
PostData(1,fOutput);
- } //this is for MC JPSI -> ee from B-hadrons
- } //this is for MC JPSI -> ee
- } //this is for MC displaced JPSI
- } //this is for MC JPSI
- } // end of check if JPSI is signal
+ } //this is for MC JPSI -> ee from B-hadrons
+ } //this is for MC JPSI -> ee
+ }
+ } // end of check if JPSI is signal
}
}
void SetPtCuts(const Double_t ptcuts[2]);\r
void SetAODMCInfo(Bool_t OkMCInfo) { fOkAODMC = OkMCInfo;}\r
void SetMinimize(Bool_t OkMinimize) { fOkMinimize = OkMinimize;}\r
- void ReadAODMCInfo(const AliAODEvent* aodEv, const TClonesArray* inArray);\r
+ void ReadAODMCInfo(AliAODEvent* aodEv, const TClonesArray* inArray);\r
void SetPathMCfile (TString mcfilename="CsiMCfromKine_10PtBins.root") {fNameMCfile = mcfilename;}\r
TH1F* OpenLocalMCFile(TString dataDir, Int_t nbin); \r
\r
Double_t fCuts[9]; // cuts for N-tuple values selection\r
Double_t fPtCuts[2]; // Max and min pt of the candidates\r
\r
- ClassDef(AliAnalysisTaskSEBtoJPSItoEle,0); // AliAnalysisTaskSE for the reconstruction of heavy-flavour decay candidates\r
+ ClassDef(AliAnalysisTaskSEBtoJPSItoEle,1); // AliAnalysisTaskSE for the reconstruction of heavy-flavour decay candidates\r
};\r
#endif\r
fFPlus(0.),
fFMinus(0.),
fFSym(0.),
-fIntegral(0.),
+fIntegral(1.),
+fintxFunB(1.),
+fintxDecayTimeBkgPos(1.),
+fintxDecayTimeBkgNeg(1.),
+fintxDecayTimeBkgSym(1.),
+fintmMassSig(1.),
+fintxRes(1.),
+fintmMassBkg(1.),
fhCsiMC(0x0),
fMassWndHigh(0.),
fMassWndLow(0.),
//
// constructor
//
- SetCrystalBallParam(kFALSE);
+ SetCrystalBallFunction(kFALSE);
SetMassWndHigh(0.2);
SetMassWndLow(0.5);
- for(Int_t iPar = 0; iPar < 13; iPar++) fParameters[iPar] = 0.;
+ for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = 0.;
fParameters[9] = 1.;fParameters[11] = 1.;fParameters[12] = 1.;
- for(Int_t index=0; index<7; index++) fResolutionConstants[index] = 0.;
+ for(Int_t index=0; index<6; index++) fResolutionConstants[index] = 0.;
AliInfo("Instance of AliBtoJPSItoEleCDFfitFCN-class created");
}
//_________________________________________________________________________________________________
fFMinus(source.fFMinus),
fFSym(source.fFSym),
fIntegral(source.fIntegral),
+fintxFunB(source.fintxFunB),
+fintxDecayTimeBkgPos(source.fintxDecayTimeBkgPos),
+fintxDecayTimeBkgNeg(source.fintxDecayTimeBkgNeg),
+fintxDecayTimeBkgSym(source.fintxDecayTimeBkgSym),
+fintmMassSig(source.fintmMassSig),
+fintxRes(source.fintxRes),
+fintmMassBkg(source.fintmMassBkg),
fhCsiMC(source.fhCsiMC),
fMassWndHigh(source.fMassWndHigh),
fMassWndLow(source.fMassWndLow),
//
// Copy constructor
//
- for(Int_t iPar = 0; iPar < 13; iPar++) fParameters[iPar] = source.fParameters[iPar];
- for(Int_t index=0; index<7; index++) fResolutionConstants[index] = source.fResolutionConstants[index];
+ for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = source.fParameters[iPar];
+ for(Int_t index=0; index<6; index++) fResolutionConstants[index] = source.fResolutionConstants[index];
}
//_________________________________________________________________________________________________
AliBtoJPSItoEleCDFfitFCN& AliBtoJPSItoEleCDFfitFCN::operator=(const AliBtoJPSItoEleCDFfitFCN& source)
fFMinus = source.fFMinus;
fFSym = source.fFSym;
fIntegral = source.fIntegral;
+ fintxFunB = source.fintxFunB;
+ fintxDecayTimeBkgPos = source.fintxDecayTimeBkgPos;
+ fintxDecayTimeBkgNeg = source.fintxDecayTimeBkgNeg;
+ fintxDecayTimeBkgSym = source.fintxDecayTimeBkgSym;
+ fintmMassSig = source.fintmMassSig;
+ fintxRes = source.fintxRes;
+ fintmMassBkg = source.fintmMassBkg;
fhCsiMC = source.fhCsiMC;
fCrystalBallParam = source.fCrystalBallParam;
- for(Int_t iPar = 0; iPar < 13; iPar++) fParameters[iPar] = source.fParameters[iPar];
- for(Int_t index=0; index<7; index++) fResolutionConstants[index] = source.fResolutionConstants[index];
+ for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = source.fParameters[iPar];
+ for(Int_t index=0; index<6; index++) fResolutionConstants[index] = source.fResolutionConstants[index];
return *this;
}
AliBtoJPSItoEleCDFfitFCN::~AliBtoJPSItoEleCDFfitFCN()
{
//
- // Default destructor
+ // Default destructor
//
delete fhCsiMC;
- for(Int_t iPar = 0; iPar < 13; iPar++) fParameters[iPar] = 0.;
- for(Int_t index=0; index<7; index++) fResolutionConstants[index] = 0.;
+ for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = 0.;
+ for(Int_t index=0; index<6; index++) fResolutionConstants[index] = 0.;
}
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t* pseudoproperdecaytime,
const Double_t* invariantmass, const Int_t ncand)
{
-//
-// This function evaluates the Likelihood fnction
-// It returns the -Log(of the likelihood function)
-//
+ //
+ // This function evaluates the Likelihood fnction
+ // It returns the -Log(of the likelihood function)
+ //
Double_t f = 0.;
Double_t ret = 0.;
for(Int_t i=0; i < ncand; i++) {
f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i]);
if(f < 0.) {
- //AliWarning("One negative contributors in the Log(Likely) ! ");
continue;
}
ret+=-1.*TMath::Log(f);
//
// Sets array of FCN parameters
//
- for(Int_t index = 0; index < 13; index++) fParameters[index] = parameters[index];
+ for(Int_t index = 0; index < 16; index++) fParameters[index] = parameters[index];
}
//_________________________________________________________________________________________________
void AliBtoJPSItoEleCDFfitFCN::ComputeIntegral()
{
-//
-// this function compute the integral of the likelihood function
-// (theoretical function) in order to normalize it to unity
-//
- Double_t np = 50.0; //number of integration steps
- Double_t stepm;Double_t stepx; //integration step width in variable m,x
- Double_t mx;Double_t xx;
+ //
+ // this function compute the integral of the likelihood function
+ // (theoretical function) in order to normalize it to unity
+ //
+ Double_t np = 100.0; //number of integration steps
+ Double_t npres = 200.0; //number of integration steps for the resolution function
+ Double_t npm = 200.;
+ Double_t stepm;Double_t stepx;Double_t stepxres; //integration step width in variable m,x
+ Double_t mx=0.;Double_t xprime=0.;
Double_t xlow = -4000.; Double_t xup = 4000.;
Double_t i; Double_t j;
Double_t sumx = 0.0;Double_t intx = 0.0;Double_t intm = 0.0;
- stepm = (fMassWndHigh-fMassWndLow)/np;
+ stepm = (fMassWndHigh-fMassWndLow)/npm;
stepx = (xup-xlow)/np;
-
- for(i = 1.0; i <= np; i++) {
- Double_t summ = 0.0;
- xx = xlow + (i - .5)*stepx;
- for(j = 1.0; j <= np/2; j++) {
- mx = fMassWndLow + (j - .5)*stepm;
- summ += EvaluateCDFfunc(xx,mx);
- mx = fMassWndHigh - (j - .5)*stepm;
- summ += EvaluateCDFfunc(xx,mx);
- }
- intm = summ*stepm;
- sumx += intm;
- }
- intx = sumx*stepx;
- SetIntegral(intx);
+ stepxres = (xup-xlow)/npres;
+
+// compute integrals for all the terms
+
+ Double_t iRes;
+ Double_t intxRes = 0.0;
+ Double_t sumxRes = 0.0;
+ for(iRes = 1.0; iRes <= npres/2.; iRes++) {
+ xprime = xlow + (iRes - .5)*stepxres;
+ sumxRes += ResolutionFunc(xprime);
+ xprime = xup - (iRes - .5)*stepxres;
+ sumxRes += ResolutionFunc(xprime);
+ }
+ intxRes = sumxRes*stepxres;
+ SetIntegralRes(intxRes);
+
+//
+ Double_t iFunB;
+ Double_t intxFunB = 0.0;
+ Double_t sumxFunB = 0.0;
+ for(iFunB = 1.0; iFunB <= np/2; iFunB++) {
+ xprime = xlow + (iFunB - .5)*stepx;
+ sumxFunB += FunB(xprime);
+ xprime = xup - (iFunB - .5)*stepx;
+ sumxFunB += FunB(xprime);
+ }
+ intxFunB = sumxFunB*stepx;
+ SetIntegralFunB(intxFunB);
+
+//
+ Double_t iDecayTimeBkgPos;
+ Double_t intxDecayTimeBkgPos = 0.0;
+ Double_t sumxDecayTimeBkgPos = 0.0;
+ for(iDecayTimeBkgPos = 1.0; iDecayTimeBkgPos <= np/2; iDecayTimeBkgPos++) {
+ xprime = xlow + (iDecayTimeBkgPos - .5)*stepx;
+ sumxDecayTimeBkgPos += FunBkgPos(xprime);
+ xprime = xup - (iDecayTimeBkgPos - .5)*stepx;
+ sumxDecayTimeBkgPos += FunBkgPos(xprime);
+ }
+ intxDecayTimeBkgPos = sumxDecayTimeBkgPos*stepx;
+ SetIntegralBkgPos(intxDecayTimeBkgPos);
+
+//
+ Double_t iDecayTimeBkgNeg;
+ Double_t intxDecayTimeBkgNeg = 0.0;
+ Double_t sumxDecayTimeBkgNeg = 0.0;
+ for(iDecayTimeBkgNeg = 1.0; iDecayTimeBkgNeg<= np/2; iDecayTimeBkgNeg++) {
+ xprime = xlow + (iDecayTimeBkgNeg - .5)*stepx;
+ sumxDecayTimeBkgNeg += FunBkgNeg(xprime);
+ xprime = xup - (iDecayTimeBkgNeg - .5)*stepx;
+ sumxDecayTimeBkgNeg += FunBkgNeg(xprime);
+ }
+ intxDecayTimeBkgNeg = sumxDecayTimeBkgNeg*stepx;
+ SetIntegralBkgNeg(intxDecayTimeBkgNeg);
+//
+ Double_t iDecayTimeBkgSym;
+ Double_t intxDecayTimeBkgSym = 0.0;
+ Double_t sumxDecayTimeBkgSym = 0.0;
+ for(iDecayTimeBkgSym = 1.0; intxDecayTimeBkgSym <= np/2; intxDecayTimeBkgSym++) {
+ xprime = xlow + (intxDecayTimeBkgSym - .5)*stepx;
+ sumxDecayTimeBkgSym += FunBkgSym(xprime);
+ xprime = xup - (intxDecayTimeBkgSym - .5)*stepx;
+ sumxDecayTimeBkgSym += FunBkgSym(xprime);
+ }
+ intxDecayTimeBkgSym = sumxDecayTimeBkgSym*stepx;
+ SetIntegralBkgSym(intxDecayTimeBkgSym);
+
+//
+ Double_t iMassSig;
+ Double_t intmMassSig = 0.0;
+ Double_t summMassSig = 0.0;
+ for(iMassSig = 1.0; iMassSig<= npm/2.; iMassSig++) {
+ mx = fMassWndLow + (iMassSig - .5)*stepm;
+ summMassSig += EvaluateCDFInvMassSigDistr(mx);
+ mx = fMassWndHigh - (iMassSig - .5)*stepm;
+ summMassSig += EvaluateCDFInvMassSigDistr(mx);
+ }
+ intmMassSig = summMassSig*stepm;
+ SetIntegralMassSig(intmMassSig);
+
+//
+ Double_t iMassBkg;
+ Double_t intmMassBkg = 0.0;
+ Double_t summMassBkg = 0.0;
+ for(iMassBkg = 1.0; iMassBkg <= npm/2.; iMassBkg++) {
+ mx = fMassWndLow + (iMassBkg - .5)*stepm;
+ summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
+ mx = fMassWndHigh - (iMassBkg - .5)*stepm;
+ summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
+ }
+ intmMassBkg = summMassBkg*stepm;
+ SetIntegralMassBkg(intmMassBkg);
+
+//
+// Compute integral of the whole distribution function
+//
+ for(i = 1.0; i <= np; i++) {
+ Double_t summ = 0.0;
+ xprime = xlow + (i - .5)*stepx;
+ for(j = 1.0; j <= npm/2; j++) {
+ mx = fMassWndLow + (j - .5)*stepm;
+ summ += EvaluateCDFfunc(xprime,mx);
+ mx = fMassWndHigh - (j - .5)*stepm;
+ summ += EvaluateCDFfunc(xprime,mx);
+ }
+ intm = summ*stepm;
+ sumx += intm;
+ }
+ intx = sumx*stepx;
+ SetIntegral(intx);
+
}
//_________________________________________________________________________________________________
void AliBtoJPSItoEleCDFfitFCN::PrintStatus()
{
-//
-// Print the parameters of the fits
-//
+ //
+ // Print the parameters of the fits
+ //
printf("\n");
printf("actual value of fRadius---------------------------------------->> | %f \n", GetRadius());
printf("actual value of fTheta ---------------------------------------->> | %f \n", GetTheta());
printf("actual value of fCrystalBallNexp ------------------------------>> | %f \n", GetCrystalBallNexp());
printf("actual value of fCrystalBallSigma ----------------------------->> | %f \n", GetCrystalBallSigma());
printf("actual value of fCrystalBallAlpha ----------------------------->> | %f \n", GetCrystalBallAlpha());
+ printf("actual value of fCrystalBallNorm ----------------------------->> | %f \n", GetCrystalBallNorm());
}else{
printf("actual value of fMpv ------------------------------------------>> | %f \n", GetCrystalBallMmean());
printf("actual value of fConstRovL ------------------------------------>> | %f \n", GetCrystalBallNexp());
printf("actual value of fSigmaL --------------------------------------->> | %f \n", GetCrystalBallSigma());
printf("actual value of fSigmaR --------------------------------------->> | %f \n", GetCrystalBallAlpha());
}
+ printf("actual value of fSigmaResol ----------------------------------->> | %f \n", GetSigmaResol());
+ printf("actual value of fNResol --------------------------------------->> | %f \n", GetNResol());
printf("\n");
- printf("Actual value of normalization integral for FCN ---------------->> | %f \n", GetIntegral());
+ printf("Actual value of normalization integral for FunB ------------------->> | %f \n", GetIntegralFunB());
+ printf("Actual value of normalization integral for BkgPos ----------------->> | %f \n", GetIntegralBkgPos());
+ printf("Actual value of normalization integral for BkgNeg ----------------->> | %f \n", GetIntegralBkgNeg());
+ printf("Actual value of normalization integral for BkgSym ----------------->> | %f \n", GetIntegralBkgSym());
+ printf("Actual value of normalization integral for MassSig ---------------->> | %f \n", GetIntegralMassSig());
+ printf("Actual value of normalization integral for MassBkg ---------------->> | %f \n", GetIntegralMassBkg());
+ printf("Actual value of normalization integral for Resolution ------------->> | %f \n", GetIntegralRes());
+ printf("Actual value of normalization integral for FCN -------------------->> | %f \n", GetIntegral());
+
printf("\n");
}
//_________________________________________________________________________________________________
-void AliBtoJPSItoEleCDFfitFCN::SetResolutionConstants(Int_t BinNum)
+void AliBtoJPSItoEleCDFfitFCN::SetResolutionConstants()
{
-// This method must be update.
-// For the time beeing the values are hard-wired.
-// Implementations have to be done to set the values from outside (e.g. from a ConfigHF file)
-//
- switch(BinNum){
-
- case(kallpt):
- fResolutionConstants[0] = 0.326 ; fResolutionConstants[2] = 0.3622*535.9 ; fResolutionConstants[4] = 1.171*5535.9 ;//from fit integrated in pt
- fResolutionConstants[1] = 0.0998*535.9; fResolutionConstants[3] = 0.1072 ; fResolutionConstants[5] = 0.04115 ; fResolutionConstants[6] = 1e-04;
- break;
- case(kptbin1):
- fResolutionConstants[0] = 0.326 ; fResolutionConstants[2] = 0.3622*1087 ; fResolutionConstants[4] = 1.171*1087 ;//from fit integrated in pt
- fResolutionConstants[1] = 0.04253*1087 ; fResolutionConstants[3] = 0.1482 ; fResolutionConstants[5] = 0.09778 ; fResolutionConstants[6] = 3.773e-04;
- break;
- case(kptbin2):
- fResolutionConstants[0] = 0.326 ; fResolutionConstants[2] = 0.3622*661.5 ; fResolutionConstants[4] = 1.171*661.5 ;//from fit integrated in pt
- fResolutionConstants[1] = 0.1*661.5 ; fResolutionConstants[3] = 0.2809 ; fResolutionConstants[5] = 0.09771 ; fResolutionConstants[6] = 1.916e-04;
- break;
- case(kptbin3):
- fResolutionConstants[0] = 0.326 ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8 ;//from fit integrated in pt
- fResolutionConstants[1] = 0.1578*502.8 ; fResolutionConstants[3] = 0.3547 ; fResolutionConstants[5] = 0.09896 ; fResolutionConstants[6] = 5.241e-04;
- break;
- case(kptbin4):
- fResolutionConstants[0] = 0.326 ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8 ;//from fit integrated in pt
- fResolutionConstants[1] = 0.2048*415.9 ; fResolutionConstants[3] = 0.4265 ; fResolutionConstants[5] = 0.09597 ; fResolutionConstants[6] = 6.469e-04;
- break;
- case(kptbin5):
- fResolutionConstants[0] = 0.326 ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8 ;//from fit integrated in pt
- fResolutionConstants[1] = 0.2219*379.7 ; fResolutionConstants[3] = 0.5414 ; fResolutionConstants[5] = 0.07506 ; fResolutionConstants[6] = 7.465e-04;
- break;
- case(kptbin6):
- fResolutionConstants[0] = 0.326 ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8 ;//from fit integrated in pt
- fResolutionConstants[1] = 0.2481*307. ; fResolutionConstants[3] = 0.8073 ; fResolutionConstants[5] = 0.09664 ; fResolutionConstants[6] = 8.481e-04;
- break;
- case(kptbin7):
- fResolutionConstants[0] = 0.326 ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8 ;//from fit integrated in pt
- fResolutionConstants[1] = 0.262*283.5 ; fResolutionConstants[3] = 0.9639 ; fResolutionConstants[5] = 0.07943 ; fResolutionConstants[6] = 6.873e-04;
- break;
- case(kptbin8):
- fResolutionConstants[0] = 0.326 ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8 ;//from fit integrated in pt
- fResolutionConstants[1] = 0.4514*204.8 ; fResolutionConstants[3] = 0.98 ; fResolutionConstants[5] = 0.1192 ; fResolutionConstants[6] = 8.646e-04;
- break;
- case(kptbin9):
- fResolutionConstants[0] = 0.326 ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8 ;//from fit integrated in pt
- fResolutionConstants[1] = 0.525*181. ; fResolutionConstants[3] = 0.99 ; fResolutionConstants[5] = 0.1097 ; fResolutionConstants[6] = 9.637e-04;
- break;
- }
+ //
+ // This method must be update:
+ // for the time beeing the values are hard-wired.
+ // Implementations have to be done to set the values from outside
+ // (e.g. from a ConfigHF file) starting from an indipendent fit
+ // of primary JPSI distribution.
+ //
+
+ fResolutionConstants[0] = 8.; // mean sigma2/sigma1
+ fResolutionConstants[1] = 0.1675; // mean Integral2/Integral1
+ fResolutionConstants[2] = 1374.; // sigma2
+ fResolutionConstants[3] = 0.001022; // N2
+ fResolutionConstants[4] = 686.6; // mu2
}
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m) const
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m) const
{
- return EvaluateCDFDecayTimeSigDistr(x)*EvaluateCDFInvMassSigDistr(m);
+ return EvaluateCDFDecayTimeSigDistr(x)*(EvaluateCDFInvMassSigDistr(m)/fintmMassSig);
}
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x) const
{
-//
-// Implementation of the Background part of the Likelyhood function
-//
-//
+ //
+ // Implementation of the Background part of the Likelyhood function
+ //
+
Double_t retvalue = 0.;
- retvalue = fParameters[7]*FunB(x) + (1. - fParameters[7])*FunP(x);
+ Double_t FunBnorm = FunB(x)/fintxFunB;
+ Double_t FunPnorm = ResolutionFunc(x)/fintxRes;
+ retvalue = fParameters[7]*FunBnorm + (1. - fParameters[7])*FunPnorm;
return retvalue;
}
//_________________________________________________________________________________________________
// Parametrization of signal part invariant mass distribution
// It can be either Crystal Ball function or sum of two Landau
//
+
Double_t fitval = 0.;
if(fCrystalBallParam){
- Double_t fitvalCB = 0.;
- Double_t normFactorCB = 1.;
- Double_t arg = (m - fParameters[9])/fParameters[11];
- Double_t numfactor = fParameters[10];
- Double_t denomfactor = numfactor - TMath::Abs(fParameters[12]) - arg;
-
- if(arg <= -1*TMath::Abs(fParameters[12])){
- Double_t exponent = fParameters[10]*TMath::Abs(fParameters[12]);
- Double_t numer = TMath::Exp(-0.5*fParameters[12]*fParameters[12])*TMath::Power(numfactor,exponent);
- Double_t denom = TMath::Power(denomfactor,exponent);
- fitvalCB += numer/denom;
- }
- if(arg > -1*TMath::Abs(fParameters[12])){
- fitvalCB += TMath::Exp(-0.5*arg*arg);
+ Double_t t = (m-fParameters[9])/fParameters[11]; ;
+ if (fParameters[12] < 0) t = -t;
+
+ Double_t absAlpha = TMath::Abs((Double_t)fParameters[12]);
+
+ if (t >= -absAlpha) {
+ return fParameters[13]*TMath::Exp(-0.5*t*t);
}
- fitval = normFactorCB*fitvalCB;
- return fitval;
+ else {
+ Double_t a = TMath::Power(fParameters[10]/absAlpha,fParameters[10])* TMath::Exp(-0.5*absAlpha*absAlpha);
+ Double_t b= fParameters[10]/absAlpha - absAlpha;
+ fitval = (fParameters[13]*a/TMath::Power(b - t, fParameters[10]));
+ return fitval;
+ }
}else{
Double_t t=-1*m;
Double_t tmpv=-1*fParameters[9];
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::FunB(Double_t x) const
{
-//
-// Parameterisation of the fit function for the x part of the Background
-//
- Double_t np = 50.0;
- Double_t sc = 100.;
+ //
+ // Parameterisation of the fit function for the x part of the Background
+ //
+
+ Double_t np = 100.0;
+ Double_t sc = 10.;
Double_t sigma3 = fResolutionConstants[2];
- Double_t xx;
+ Double_t xprime;
Double_t sum = 0.0;
Double_t xlow,xupp;
Double_t step;
xlow = x - sc * sigma3 ;
xupp = x + sc * sigma3 ;
step = (xupp-xlow) / np;
- for(i=1.0; i<=np/2; i++) {
- xx = xlow + (i-.5) * step;
- sum += CsiMC(xx) * ResolutionFunc(xx-x);
+ Double_t CsiMCxprime = 0.;
+ Double_t Resolutionxdiff = 0.;
+ Double_t xdiff = 0.;
+
+ for(i=1.0; i<=np; i++) {
+
+ xprime = xlow + (i-.5) * step;
+ CsiMCxprime = CsiMC(xprime);
+ xdiff = xprime - x;
+ Resolutionxdiff = ResolutionFunc(xdiff)/fintxRes; // normalized value
+ sum += CsiMCxprime * Resolutionxdiff;
- xx = xupp - (i-.5) * step;
- sum += CsiMC(xx) * ResolutionFunc(xx-x);
}
- return (step * sum) ;
+
+ return step * sum ;
}
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::FunP(Double_t x) const
-//
-// Parameterisation of the Prompt part for the x distribution
-//
{
- return ResolutionFunc(x);
+ //
+ // Parameterisation of the Prompt part for the x distribution
+ //
+
+ return ResolutionFunc(x)/fintxRes; // normalized value
}
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::CsiMC(Double_t x) const
{
-// Distribution (template) of the x distribution for the x variable
-// for the J/psi coming from Beauty hadrons
-//
+ //
+ // Distribution (template) of the x distribution for the x variable
+ // for the J/psi coming from Beauty hadrons
+ //
+
Double_t returnvalue = 0.;
returnvalue = fhCsiMC->GetBinContent(fhCsiMC->FindBin(x));
return returnvalue;
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m) const
{
-//
-// Return the part of the likelihood function for the background hypothesis
-//
- return EvaluateCDFDecayTimeBkgDistr(x)*EvaluateCDFInvMassBkgDistr(m);
+ //
+ // Return the part of the likelihood function for the background hypothesis
+ //
+
+ return EvaluateCDFDecayTimeBkgDistr(x)*(EvaluateCDFInvMassBkgDistr(m)/fintmMassBkg);
}
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x) const
{
-// it returns the value of the probability to have a given x for the background
-//
-//
- Double_t ret = (1 - fParameters[0]*fParameters[0])*ResolutionFunc(x);
- ret += FunBkgPos(x);
- ret += FunBkgNeg(x);
- ret += FunBkgSym(x);
+ //
+ // it returns the value of the probability to have a given x for the background
+ //
+
+ Double_t ret = (1. - TMath::Power(fParameters[0],2.))*(ResolutionFunc(x)/fintxRes)
+ + TMath::Power(fParameters[0]*TMath::Cos(fParameters[1]),2.)*
+ (FunBkgPos(x)/fintxDecayTimeBkgPos)
+ + TMath::Power(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]),2.)*
+ (FunBkgNeg(x)/fintxDecayTimeBkgNeg)
+ + TMath::Power(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]),2.)*
+ (FunBkgSym(x)/fintxDecayTimeBkgSym);
return ret;
}
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistr(Double_t m) const
{
-//
-// it returns the value of the probability to have a given mass for the background
-//
- return 1/(fMassWndHigh-fMassWndLow)+fParameters[6]*(m-(fMassWndHigh+fMassWndLow)/2);
+ //
+ // it returns the value of the probability to have a given mass for the background
+ //
+
+ return 1/(fMassWndHigh-fMassWndLow) +
+ fParameters[6] * m -
+ fParameters[6] * ((fMassWndHigh+fMassWndLow)/2);
}
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x) const
{
-//
-// exponential with positive slopes for the background part (x)
-//
- Double_t np = 50.0;
- Double_t sc = 100.;
+ //
+ // exponential with positive slopes for the background part (x)
+ //
+
+ Double_t np = 100.0;
+ Double_t sc = 10.;
Double_t sigma3 = fResolutionConstants[2];
- Double_t xx;
+ Double_t xprime;
Double_t sum = 0.0;
Double_t xlow,xupp;
Double_t step;
xlow = x - sc * sigma3 ;
xupp = x + sc * sigma3 ;
step = (xupp-xlow) / np;
+
for(i=1.0; i<=np/2; i++) {
- xx = xlow + (i-.5) * step;
- if (xx > 0) sum += (fParameters[0]*TMath::Cos(fParameters[1]))*(fParameters[0]*TMath::Cos(fParameters[1]))*fParameters[3]*TMath::Exp(-1*xx*fParameters[3]) * ResolutionFunc(xx-x);
+
+ xprime = xlow + (i-.5) * step;
+ if (xprime > 0) {sum += fParameters[3] * TMath::Exp(-1*xprime*fParameters[3]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
- xx = xupp - (i-.5) * step;
- if (xx > 0) sum += fParameters[0]*TMath::Cos(fParameters[1])*(fParameters[0]*TMath::Cos(fParameters[1]))*fParameters[3]*TMath::Exp(-1*xx*fParameters[3]) * ResolutionFunc(xx-x);
+ xprime = xupp - (i-.5) * step;
+ if (xprime > 0) {sum += fParameters[3] * TMath::Exp(-1*xprime*fParameters[3]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
+
}
- return (step * sum) ;
+
+ return step * sum ;
}
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x) const
{
-//
-// exponential with negative slopes for the background part (x)
-//
- Double_t np = 50.0;
- Double_t sc = 100.;
+ //
+ // exponential with negative slopes for the background part (x)
+ //
+
+ Double_t np = 100.0;
+ Double_t sc = 10.;
Double_t sigma3 = fResolutionConstants[2];
- Double_t xx;
+ Double_t xprime;
Double_t sum = 0.0;
Double_t xlow,xupp;
Double_t step;
xlow = x - sc * sigma3 ;
xupp = x + sc * sigma3 ;
step = (xupp-xlow) / np;
+
for(i=1.0; i<=np/2; i++) {
- xx = xlow + (i-.5) * step;
- if (xx < 0) sum += (fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]))*
- fParameters[4]*TMath::Exp(xx*fParameters[4]) * ResolutionFunc(xx-x);
- xx = xupp - (i-.5) * step;
- if (xx < 0) sum += (fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]))*
- fParameters[4]*TMath::Exp(xx*fParameters[4]) * ResolutionFunc(xx-x);
+
+ xprime = xlow + (i-.5) * step;
+ if (xprime < 0) {sum += fParameters[4] * TMath::Exp(xprime*fParameters[4]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
+
+ xprime = xupp - (i-.5) * step;
+ if (xprime < 0) {sum += fParameters[4] * TMath::Exp(xprime*fParameters[4]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
}
- return (step * sum) ;
+
+ return step * sum ;
}
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x) const
{
-//
-// exponential with both positive and negative slopes for the background part (x)
-//
- Double_t np = 50.0;
- Double_t sc = 100.;
+ //
+ // exponential with both positive and negative slopes for the background part (x)
+ //
+
+ Double_t np = 100.0;
+ Double_t sc = 10.;
Double_t sigma3 = fResolutionConstants[2];
- Double_t xx;
+ Double_t xprime;
Double_t sum1 = 0.0;
Double_t sum2 = 0.0;
Double_t xlow,xupp;
xlow = x - sc * sigma3 ;
xupp = x + sc * sigma3 ;
step = (xupp-xlow) / np;
+
for(i=1.0; i<=np/2; i++) {
- xx = xlow + (i-.5) * step;
- if (xx > 0) sum1 += 0.5*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*
- fParameters[5]*TMath::Exp(-1*xx*fParameters[5]) * ResolutionFunc(xx-x);
+
+ xprime = xlow + (i-.5) * step;
+ if (xprime > 0) {sum1 += 0.5 * fParameters[5]*TMath::Exp(-1*xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum1 = 0.;}
- xx = xupp - (i-.5) * step;
- if (xx > 0) sum1 += 0.5*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*
- fParameters[5]*TMath::Exp(-1*xx*fParameters[5]) * ResolutionFunc(xx-x);
+ xprime = xupp - (i-.5) * step;
+ if (xprime > 0) {sum1 += 0.5 * fParameters[5]*TMath::Exp(-1*xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum1 = 0.;}
}
+
for(i=1.0; i<=np/2; i++) {
- xx = xlow + (i-.5) * step;
- if (xx < 0) sum2 += 0.5*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*
- fParameters[5]*TMath::Exp(xx*fParameters[5]) * ResolutionFunc(xx-x);
- xx = xupp - (i-.5) * step;
- if (xx < 0) sum2 += 0.5*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*
- fParameters[5]*TMath::Exp(xx*fParameters[5]) * ResolutionFunc(xx-x);
+ xprime = xlow + (i-.5) * step;
+ if (xprime < 0) {sum2 += 0.5 * fParameters[5]*TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum2 = 0.;}
+
+ xprime = xupp - (i-.5) * step;
+ if (xprime < 0) {sum2 += 0.5 * fParameters[5]*TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum2 = 0.;}
}
+
return step*(sum1 + sum2) ;
}
//_________________________________________________________________________________________________
Double_t AliBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x) const
{
//
- //parametrization with 2 gaus + 1 exponential + 1 constant
- //
- Double_t arg=0;
- arg=x/fResolutionConstants[1];
- Double_t ret=TMath::Exp(-0.5*arg*arg);
- arg=x/fResolutionConstants[2];
- ret+=fResolutionConstants[3]*TMath::Exp(-0.5*arg*arg);
- arg=x/fResolutionConstants[4];
- if(x > 0) { ret+=fResolutionConstants[5]*TMath::Exp(-arg); }
- if(x <= 0) { ret+=fResolutionConstants[5]*TMath::Exp(arg); }
- return fResolutionConstants[0]*(ret + fResolutionConstants[6]);
+ // parametrization with 2 gaus
+ //
+
+ Double_t ret = 0.;
+ Double_t x1 = x;
+ Double_t x2 = x;
+ //Double_t mean1 = 0.;
+ Double_t mean2 = fResolutionConstants[4];
+ Double_t sigma1 = fParameters[14];
+ Double_t sigma2 = fResolutionConstants[2];
+ Double_t n1 = fParameters[15];
+ Double_t n2 = fResolutionConstants[3];
+ Double_t arg1 = x1/sigma1;
+ Double_t arg2 = (x2-mean2)/sigma2;
+ Double_t sqrt2Pi = TMath::Sqrt(2*TMath::Pi());
+
+ ret = n2*((n1/n2)*TMath::Exp(-0.5*arg1*arg1) + TMath::Exp(-0.5*arg2*arg2));
+
+ return ret/(sqrt2Pi*(n1*sigma1+n2*sigma2));//return value is normalized
+
}
+
#include <TDatabasePDG.h>\r
#include "TH1F.h"\r
#include "TMath.h"\r
+#include "TRandom3.h"\r
\r
enum IntegralType {kBkg, \r
kBkgNorm, \r
Double_t EvaluateLikelihood(const Double_t* pseudoproperdecaytime,\r
const Double_t* invariantmass, const Int_t ncand);\r
\r
- Double_t GetFPlus() const {return fFPlus;}\r
- Double_t GetFMinus() const {return fFMinus;}\r
- Double_t GetFSym() const {return fFSym;}\r
- Double_t GetRadius() const {return fParameters[0];}\r
- Double_t GetTheta() const {return fParameters[1];}\r
- Double_t GetPhi() const {return fParameters[2];}\r
- Double_t GetLamPlus() const {return fParameters[3];}\r
- Double_t GetLamMinus() const {return fParameters[4];}\r
- Double_t GetLamSym() const {return fParameters[5];}\r
- Double_t GetMassSlope() const {return fParameters[6];}\r
- Double_t GetFractionJpsiFromBeauty() const {return fParameters[7];}\r
- Double_t GetFsig() const {return fParameters[8];}\r
- Double_t GetCrystalBallMmean() const {return fParameters[9];}\r
- Double_t GetCrystalBallNexp() const {return fParameters[10];}\r
- Double_t GetCrystalBallSigma() const {return fParameters[11];}\r
- Double_t GetCrystalBallAlpha() const {return fParameters[12];}\r
- Double_t GetIntegral() const {return fIntegral;}\r
- Bool_t GetCrystalBallParam() const {return fCrystalBallParam;}\r
+ Double_t GetFPlus() const { return fFPlus; }\r
+ Double_t GetFMinus() const { return fFMinus; }\r
+ Double_t GetFSym() const { return fFSym; }\r
+ Double_t GetRadius() const { return fParameters[0]; }\r
+ Double_t GetTheta() const { return fParameters[1]; }\r
+ Double_t GetPhi() const { return fParameters[2]; }\r
+ Double_t GetLamPlus() const { return fParameters[3]; }\r
+ Double_t GetLamMinus() const { return fParameters[4]; }\r
+ Double_t GetLamSym() const { return fParameters[5]; }\r
+ Double_t GetMassSlope() const { return fParameters[6]; }\r
+ Double_t GetFractionJpsiFromBeauty() const { return fParameters[7]; }\r
+ Double_t GetFsig() const { return fParameters[8]; }\r
+ Double_t GetCrystalBallMmean() const { return fParameters[9]; }\r
+ Double_t GetCrystalBallNexp() const { return fParameters[10]; }\r
+ Double_t GetCrystalBallSigma() const { return fParameters[11]; }\r
+ Double_t GetCrystalBallAlpha() const { return fParameters[12]; }\r
+ Double_t GetCrystalBallNorm() const { return fParameters[13]; }\r
+ Double_t GetSigmaResol() const { return fParameters[14]; }\r
+ Double_t GetNResol() const { return fParameters[15]; }\r
+ Double_t GetIntegral() const { return fIntegral; }\r
+ Double_t GetIntegralFunB() const { return fintxFunB; }\r
+ Double_t GetIntegralBkgPos() const { return fintxDecayTimeBkgPos; }\r
+ Double_t GetIntegralBkgNeg() const { return fintxDecayTimeBkgNeg; }\r
+ Double_t GetIntegralBkgSym() const { return fintxDecayTimeBkgSym; }\r
+ Double_t GetIntegralMassSig() const { return fintmMassSig; }\r
+ Double_t GetIntegralRes() const { return fintxRes; }\r
+ Double_t GetIntegralMassBkg() const { return fintmMassBkg; }\r
+ Bool_t GetCrystalBallParam() const { return fCrystalBallParam; }\r
\r
void SetFPlus(Double_t plus) {fFPlus = plus;}\r
void SetFMinus(Double_t minus) {fFMinus = minus;}\r
void SetCrystalBallNexp(Double_t CrystalBallNexp) {fParameters[10] = CrystalBallNexp;}\r
void SetCrystalBallSigma(Double_t CrystalBallSigma) {fParameters[11] = CrystalBallSigma;}\r
void SetCrystalBallAlpha(Double_t CrystalBallAlpha) {fParameters[12] = CrystalBallAlpha;}\r
+ void SetCrystalBallNorm(Double_t CrystalBallNorm) {fParameters[13] = CrystalBallNorm;}\r
+ void SetSigmaResol(Double_t SigmaResol) {fParameters[14] = SigmaResol;}\r
+ void SetNResol(Double_t NResol) {fParameters[15] = NResol;}\r
\r
void SetAllParameters(const Double_t* parameters);\r
- void SetIntegral(Double_t integral) {fIntegral = integral;}\r
+ void SetIntegral(Double_t integral) { fIntegral = integral; }\r
+ void SetIntegralFunB(Double_t integral) { fintxFunB = integral; }\r
+ void SetIntegralBkgPos(Double_t integral) { fintxDecayTimeBkgPos = integral; }\r
+ void SetIntegralBkgNeg(Double_t integral) { fintxDecayTimeBkgNeg = integral; }\r
+ void SetIntegralBkgSym(Double_t integral) { fintxDecayTimeBkgSym = integral; }\r
+ void SetIntegralMassSig(Double_t integral) { fintmMassSig = integral; }\r
+ void SetIntegralRes(Double_t integral) { fintxRes = integral; }\r
+ void SetIntegralMassBkg(Double_t integral) { fintmMassBkg = integral; }\r
+\r
void SetCsiMC(const TH1F* MCtemplate) {fhCsiMC = (TH1F*)MCtemplate->Clone("fhCsiMC");}\r
- void SetResolutionConstants(Int_t BinNum);\r
- void SetMassWndHigh(Double_t limit) { fMassWndHigh = TDatabasePDG::Instance()->GetParticle(443)->Mass() + limit ;}//here use pdg code instead\r
- void SetMassWndLow(Double_t limit) { fMassWndLow = TDatabasePDG::Instance()->GetParticle(443)->Mass() - limit ;}//here use pdg code instead\r
- void SetCrystalBallParam(Bool_t okCB) {fCrystalBallParam = okCB;}\r
+\r
+ void SetResolutionConstants();\r
+ void SetMassWndHigh(Double_t limit) { fMassWndHigh = TDatabasePDG::Instance()->GetParticle(443)->Mass() + limit ;}\r
+ void SetMassWndLow(Double_t limit) { fMassWndLow = TDatabasePDG::Instance()->GetParticle(443)->Mass() - limit ;}\r
+ void SetCrystalBallFunction(Bool_t okCB) {fCrystalBallParam = okCB;}\r
\r
void ConvertFromSpherical() { fFPlus = TMath::Power((fParameters[0]*TMath::Cos(fParameters[1])),2.);\r
fFMinus = TMath::Power((fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2])),2.);\r
\r
private: \r
//\r
- Double_t fParameters[13]; /* par[0] = fRadius; \r
+ Double_t fParameters[16]; /* par[0] = fRadius; \r
par[1] = fTheta;\r
par[2] = fPhi;\r
par[3] = fOneOvLamPlus;\r
par[9] = fCrystalBallMmean;\r
par[10] = fCrystalBallNexp;\r
par[11] = fCrystalBallSigma;\r
- par[12] = fCrystalBallAlpha;*/\r
-\r
- Double_t fFPlus; // parameters of the log-likelihood function\r
- Double_t fFMinus; // Slopes of the x distributions of the background \r
- Double_t fFSym; // functions \r
-\r
- Double_t fIntegral; // integral values of log-likelihood function terms\r
- TH1F *fhCsiMC; // X distribution used as MC template for JPSI from B\r
- Double_t fResolutionConstants[7]; // constants for the parametrized resolution function R(X)\r
- Double_t fMassWndHigh; // JPSI Mass window higher limit\r
- Double_t fMassWndLow; // JPSI Mass window lower limit\r
- Bool_t fCrystalBallParam; // Boolean to switch to Crystall Ball parameterisation\r
+ par[12] = fCrystalBallAlpha;\r
+ par[13] = fCrystalBallNorm;\r
+ par[14] = fSigmaResol;\r
+ par[15] = fNResol; */\r
+\r
+\r
+ Double_t fFPlus; // parameters of the log-likelihood function\r
+ Double_t fFMinus; // Slopes of the x distributions of the background \r
+ Double_t fFSym; // functions \r
+\r
+ Double_t fIntegral; // integral values of log-likelihood function terms\r
+ Double_t fintxFunB;\r
+ Double_t fintxDecayTimeBkgPos;\r
+ Double_t fintxDecayTimeBkgNeg;\r
+ Double_t fintxDecayTimeBkgSym;\r
+ Double_t fintmMassSig;\r
+ Double_t fintxRes;\r
+ Double_t fintmMassBkg;\r
+ TH1F *fhCsiMC; // X distribution used as MC template for JPSI from B\r
+ Double_t fResolutionConstants[5]; // constants for the parametrized resolution function R(X)\r
+ Double_t fMassWndHigh; // JPSI Mass window higher limit\r
+ Double_t fMassWndLow; // JPSI Mass window lower limit\r
+ Bool_t fCrystalBallParam; // Boolean to switch to Crystall Ball parameterisation\r
\r
////\r
\r
\r
Double_t ResolutionFunc(Double_t x) const ;\r
//\r
- ClassDef (AliBtoJPSItoEleCDFfitFCN,0); // Unbinned log-likelihood fit \r
+ ClassDef (AliBtoJPSItoEleCDFfitFCN,1); // Unbinned log-likelihood fit \r
\r
};\r
\r
//_________________________________________________________________________________________________
AliBtoJPSItoEleCDFfitHandler::AliBtoJPSItoEleCDFfitHandler():
+fIsParamFixed(16),
+fPrintStatus(kFALSE),
fUp(0),
fX(0x0),
fM(0x0),
//_________________________________________________________________________________________________
AliBtoJPSItoEleCDFfitHandler::AliBtoJPSItoEleCDFfitHandler(Double_t* decaytime,
Double_t* invariantmass, Int_t ncand) :
+fIsParamFixed(16),
+fPrintStatus(kFALSE),
fUp(0),
fX(decaytime),
fM(invariantmass),
// constructor
//
AliInfo("\n+++\n+++ Minimization object AliBtoJPSItoEleCDFfitHandler created\n+++\n");
- AliInfo("\n+++\n+++ Creating AliBtoJPSItoEleCDFfitFCN object\n+++\n");
fLikely = new AliBtoJPSItoEleCDFfitFCN();
- fLikely->SetCrystalBallParam(kFALSE); //Landau selected; otherwise Crystal Ball is selected
- SetErrorDef(1.);
+ AliInfo("\n+++\n+++ CDF fit function object AliBtoJPSItoEleCDFfitFCN created\n+++\n");
+ AliInfo("Parameter 0 ----> fRadius");
+ AliInfo("Parameter 1 ----> fTheta");
+ AliInfo("Parameter 2 ----> fPhi");
+ AliInfo("Parameter 3 ----> fOneOvLamPlus");
+ AliInfo("Parameter 4 ----> fOneOvLamMinus");
+ AliInfo("Parameter 5 ----> fOneOvLamSym");
+ AliInfo("Parameter 6 ----> fMSlope");
+ AliInfo("Parameter 7 ----> fB");
+ AliInfo("Parameter 8 ----> fFsig");
+ AliInfo("Parameter 9 ----> fMmean");
+ AliInfo("Parameter 10 ----> fNexp");
+ AliInfo("Parameter 11 ----> fSigma");
+ AliInfo("Parameter 12 ----> fAlpha");
+ AliInfo("Parameter 13 ----> fNorm");
+ AliInfo("Parameter 14 ----> fSigmaResol");
+ AliInfo("Parameter 15 ----> fNResol");
AliInfo(Form("\n+++\n+++ Number of candidates ---> %d\n+++\n ", ncand));
}
//___________________________________________________________________________
// Assignment operator
//
if (this!=&c) {
- fUp = c.fUp;
- fX = c.fX;
- fM = c.fM;
- fLikely = c.fLikely;
- fNcand = c.fNcand;
+ fIsParamFixed = c.fIsParamFixed;
+ fPrintStatus = c.fPrintStatus;
+ fUp = c.fUp;
+ fX = c.fX;
+ fM = c.fM;
+ fLikely = c.fLikely;
+ fNcand = c.fNcand;
}
return *this;
}
-//___________________________________________________________________________
+//_______________________________________________________________________________________
AliBtoJPSItoEleCDFfitHandler::AliBtoJPSItoEleCDFfitHandler(const AliBtoJPSItoEleCDFfitHandler& c) :
TNamed(c),
+fIsParamFixed(c.fIsParamFixed),
+fPrintStatus(c.fPrintStatus),
fUp(c.fUp),
fX(c.fX),
fM(c.fM),
// Copy Constructor
//
}
-//________________________________________________________________________
+//_______________________________________________________________________________________
AliBtoJPSItoEleCDFfitHandler::~AliBtoJPSItoEleCDFfitHandler()
{
//
//
delete fLikely;
}
-//_________________________________________________________________________________________________
+//_______________________________________________________________________________________
Int_t AliBtoJPSItoEleCDFfitHandler::DoMinimization()
{
//
// performs the minimization
//
- static TVirtualFitter *fitter = TVirtualFitter::Fitter(this,13);
+ static TVirtualFitter *fitter = TVirtualFitter::Fitter(this,16);
fitter->SetFCN(CDFFunction);
- Double_t startingParamValues[13] =
- /* startfPlus
- startfMinus
- startfSym
- startfOneOvLamPlus
- startfOneOvLamMinus
- startfOneOvLamSym
- startfMSlope
- startfB
- startfFsig
- startfMmean
- startfNexp
- startfSigma
- startfAlpha */
- {5.00e-01,
- TMath::Pi()/4.,
- TMath::Pi()/4.,
- 2.0964360e-03,
- 4.8309180e-03,
- 1.582530e-04,
- -1.5720e-02,
- 0.1800e+00,
- 0.7000e+00,
- 3.0910e+00,
- 1.0500e+00,
- 1.4250e-02,
- 6.758e-01};
-
- fitter->SetParameter(0,"fRadius",startingParamValues[0], 0.01, 0., 1.);
- fitter->SetParameter(1,"fTheta",startingParamValues[1], 0.001, 0., TMath::Pi()/2);
- fitter->SetParameter(2,"fPhi",startingParamValues[2], 0.001, 0., TMath::Pi()/2);
- fitter->SetParameter(3,"fOneOvLamPlus",startingParamValues[3], 0.0001, 0., 5.e-01);
- fitter->SetParameter(4,"fOneOvLamMinus",startingParamValues[4], 0.0001, 0., 5.e-01);
- fitter->SetParameter(5,"fOneOvLamSym",startingParamValues[5], 0.00001, 0., 5.e-01);
- fitter->SetParameter(6,"fMSlope",startingParamValues[6], 0.001, -1.e-00, 1.e+00);
- fitter->SetParameter(7,"fB",startingParamValues[7], 0.1, 0., 1.);
- fitter->SetParameter(8,"fFsig",startingParamValues[8], 0.1, 0., 1.);
- fitter->SetParameter(9,"fMmean",startingParamValues[9], 0.1, 0., 1.e+02);
- fitter->SetParameter(10,"fNexp",startingParamValues[10], 0.1, 0., 1.e+02);
- fitter->SetParameter(11,"fSigma",startingParamValues[11], 0.001, 0., 1.e+02);
- fitter->SetParameter(12,"fAlpha",startingParamValues[12], 0.01, 0., 1.e+02);
- fitter->FixParameter(9);
-
- Double_t arglist[2]={10000,0.5};
+
+ fitter->SetParameter(0,"fRadius",fParamStartValues[0], 1.e-06, 0., 1.);
+ fitter->SetParameter(1,"fTheta",fParamStartValues[1], 1.e-06, 0.,2*TMath::Pi());
+ fitter->SetParameter(2,"fPhi",fParamStartValues[2], 1.e-06, 0.,2*TMath::Pi());
+// fitter->SetParameter(3,"fOneOvLamPlus",fParamStartValues[3], 1.e-10, 0., 5.e+01);
+// fitter->SetParameter(4,"fOneOvLamMinus",fParamStartValues[4], 1.e-10, 0., 5.e+01);
+// fitter->SetParameter(5,"fOneOvLamSym",fParamStartValues[5], 1.e-10, 0., 5.e+01);
+ fitter->SetParameter(3,"fOneOvLamPlus",fParamStartValues[3], 1.e-10, 0.0000001, 5.e+01);
+ fitter->SetParameter(4,"fOneOvLamMinus",fParamStartValues[4], 1.e-10, 0.00000001, 5.e+01);
+ fitter->SetParameter(5,"fOneOvLamSym",fParamStartValues[5], 1.e-10, 0.00000001, 5.e+01);
+ fitter->SetParameter(6,"fMSlope",fParamStartValues[6], 1.e-04, -2.5, 2.5);
+ fitter->SetParameter(7,"fB",fParamStartValues[7], 1.e-08, 0., 1.);
+ fitter->SetParameter(8,"fFsig",fParamStartValues[8], 1.e-08, 0., 1.);
+ fitter->SetParameter(9,"fMmean",fParamStartValues[9], 1.e-08, 0., 1.e+04);
+ fitter->SetParameter(10,"fNexp",fParamStartValues[10], 1.e-08, 0., 1.e+02);
+ fitter->SetParameter(11,"fSigma",fParamStartValues[11], 1.e-08, 0., 1.e+04);
+ fitter->SetParameter(12,"fAlpha",fParamStartValues[12], 1.e-08, 0., 1.e+04);
+ fitter->SetParameter(13,"fNorm",fParamStartValues[13], 1.e-08, 0., 1.e+01);
+ fitter->SetParameter(14,"fSigmaResol",fParamStartValues[14], 1.e-08, 0., 1.e+04);
+ fitter->SetParameter(15,"fNResol",fParamStartValues[15], 1.e-08, 0., 1.e+05);
+
+ for(UInt_t indexparam = 0; indexparam < 16; indexparam++){
+ if(IsParamFixed(indexparam))fitter->FixParameter((Int_t)indexparam);
+ }
+
+ Double_t arglist[2]={10000,1.0};
Int_t iret=fitter->ExecuteCommand("MIGRAD", arglist ,2);
fitter->PrintResults(4,0);
fLikely->SetAllParameters(par);
fLikely->ConvertFromSpherical();
fLikely->ComputeIntegral();
- //printf("\n+++\n+++\n+++\n");
- //fLikely->PrintStatus();
+ if(fPrintStatus)fLikely->PrintStatus();
TStopwatch t;
t.Start();
return;
}
+//_______________________________________________________________________________________
+void AliBtoJPSItoEleCDFfitHandler::SetParamStartValues(Double_t inputparamvalues[16])
+{
+ for(Int_t index=0; index < 16; index++) fParamStartValues[index] = inputparamvalues[index];
+}
+//_______________________________________________________________________________________
+void AliBtoJPSItoEleCDFfitHandler::SetResolutionConstants()
+{
+ //
+ // Sets constants for the resolution function
+ //
+ fLikely->SetResolutionConstants();
+
+}
+//_______________________________________________________________________________________
+void AliBtoJPSItoEleCDFfitHandler::SetCrystalBallFunction(Bool_t okCB)
+{
+ //
+ // Sets the CB as the parametrization for the signal invariant mass spectrum
+ // (otherwise Landau is chosen)
+ //
+ fLikely->SetCrystalBallFunction(okCB);
+}
+//_______________________________________________________________________________________
+void AliBtoJPSItoEleCDFfitHandler::SetMassWndHigh(Double_t limit)
+{
+ //
+ // Sets upper limit for the invariant mass window (under J/PSI mass peak)
+ //
+ fLikely->SetMassWndHigh(limit);
+}
+//_______________________________________________________________________________________
+void AliBtoJPSItoEleCDFfitHandler::SetMassWndLow(Double_t limit)
+{
+ //
+ // Sets lower limit for the invariant mass window (under J/PSI mass peak)
+ //
+ fLikely->SetMassWndLow(limit);
+}
+
//____________________________________________________________________________\r
\r
#include <TNamed.h>\r
+#include <TBits.h>\r
#include "TMath.h"\r
#include "AliBtoJPSItoEleCDFfitFCN.h"\r
-//#include "AliLog.h"\r
+#include "AliLog.h"\r
+\r
+class TBits; \r
\r
class AliBtoJPSItoEleCDFfitHandler : public TNamed {\r
public:\r
AliBtoJPSItoEleCDFfitHandler(const AliBtoJPSItoEleCDFfitHandler& c);\r
AliBtoJPSItoEleCDFfitHandler(Double_t* decaytime, Double_t* invariantmass, Int_t ncand);\r
~AliBtoJPSItoEleCDFfitHandler(); \r
- Double_t Up() const { return fUp*TMath::Sqrt(fLikely->GetIntegral()); }\r
+ //Double_t Up() const { return fUp*TMath::Sqrt(fLikely->GetIntegral()); }\r
+ Double_t Up() const { return fUp; }\r
void SetErrorDef(Double_t up) {fUp = up;}\r
+ void SetPrintStatus(Bool_t okPrintStatus) { fPrintStatus = okPrintStatus; } \r
+ Bool_t GetPrintStatus() { return fPrintStatus ; }\r
+ void SetParamStartValues (Double_t*);\r
+ Double_t* GetStartParamValues() { return fParamStartValues; }\r
+ TBits GetFixedParamList() { return fIsParamFixed; }\r
+ void FixParam(UInt_t param, Bool_t value) { fIsParamFixed.SetBitNumber(param,value); }\r
+ void FixAllParam(Bool_t value) { for(UInt_t par=0;par<16;par++) fIsParamFixed.SetBitNumber(par,value); }\r
+ Bool_t IsParamFixed(UInt_t param) { return fIsParamFixed.TestBitNumber(param); }\r
+ void SetResolutionConstants();\r
+ void SetCrystalBallFunction(Bool_t okCB);\r
+ void SetMassWndHigh(Double_t limit);\r
+ void SetMassWndLow(Double_t limit);\r
\r
Double_t operator()(const Double_t* par) const ;\r
void CdfFCN(Int_t & /* npar */, Double_t * /* gin */, Double_t &f, Double_t *par, Int_t /* iflag */);\r
\r
private:\r
//\r
+ TBits fIsParamFixed; //array of bits: 0 = param free; 1 = param fixed;\r
+ Bool_t fPrintStatus; //flag to enable the prit out of the algorithm at each step\r
+ Double_t fParamStartValues[16]; //array of parameters input value\r
Double_t fUp; //error definition \r
Double_t* fX; //pseudo-proper decay time X\r
Double_t* fM; //invariant mass M\r
AliBtoJPSItoEleCDFfitFCN* fLikely; //Log likelihood function\r
Int_t fNcand; //number of candidates\r
//\r
- ClassDef(AliBtoJPSItoEleCDFfitHandler,0);\r
+ ClassDef(AliBtoJPSItoEleCDFfitHandler,1);\r
\r
}; \r
#endif\r
--- /dev/null
+void FitCDFLocal() {
+
+ ///////////////////////////////////////////////////////////////////
+ //
+ // Example macro to read local N-Tuples of JPSI
+ // and bkg candidates and perform log-likelihood
+ // minimization using the minimization handler class
+ //
+ // Origin: C. Di Giglio
+ //
+ ///////////////////////////////////////////////////////////////////
+
+ Bool_t useParFiles=kFALSE;
+ gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/LoadLibraries.C");
+ LoadLibraries(useParFiles);
+
+ TH1F* hCsiMCPithya = new TH1F();
+ TH1F* hCsiMCevtgen = new TH1F();
+
+ TString datadir=".";
+ TString useFile = datadir.Data();
+ useFile.Append("/CmpBdecayKine.root"); // file with MC templates for J/psi (<-B) X distribution (in different pT bins)
+ TFile *fCmpBdecayKine = new TFile(useFile);
+ hCsiMCPithya = (TH1F*)fCmpBdecayKine->Get("PsproperPithya46GeV"); // Pithya case: 4 - 6 GeV
+ hCsiMCPithya->Scale(1./hCsiMCPithya->Integral());
+ hCsiMCevtgen = (TH1F*)fCmpBdecayKine->Get("PsproperEvtGen46GeV"); // EvtGen case: 4 - 6 GeV
+ hCsiMCevtgen->Scale(1./hCsiMCevtgen->Integral());
+
+ Double_t* x=0x0; Double_t* m=0x0; Int_t n=0;
+ AliAnalysisBtoJPSItoEle* aBtoJPSItoEle =new AliAnalysisBtoJPSItoEle();
+ Double_t paramInputValues [16] = /*
+ paramInputValues[0] ----> fRadius;
+ paramInputValues[1] ----> fTheta;
+ paramInputValues[2] ----> fPhi;
+ paramInputValues[3] ----> fOneOvLamPlus;
+ paramInputValues[4] ----> fOneOvLamMinus;
+ paramInputValues[5] ----> fOneOvLamSym;
+ paramInputValues[6] ----> fMSlope;
+ paramInputValues[7] ----> fB;
+ paramInputValues[8] ----> fFsig;
+ paramInputValues[9] ----> fMmean;
+ paramInputValues[10] ----> fNexp;
+ paramInputValues[11] ----> fSigma;
+ paramInputValues[12] ----> fAlpha;
+ paramInputValues[13] ----> fNorm;
+ paramInputValues[14] ----> fSigmaResol;
+ paramInputValues[15] ----> fNResol;
+ */
+ { 0.5,TMath::Pi()/4.,
+ TMath::Pi()/4.,
+ 0.00210,
+ 0.00480,
+ 0.,
+ 0.,
+ 0.227272,
+ 0.3981,
+ 3.09568,
+ 0.9671,
+ 0.0239,
+ 0.6599,
+ 0.04587,
+ 55.4,
+ 0.033 }; // Starting values for parameters
+
+ TFile *f= new TFile("CdfFit_OneYear.root"); // file with N-tuples for one year collected statistics (in different pT bins)
+ //TNtuple *nt=(TNtuple*)f->Get("fNtupleJPSI12GeV");
+ //TNtuple *nt=(TNtuple*)f->Get("fNtupleJPSI23GeV");
+ //TNtuple *nt=(TNtuple*)f->Get("fNtupleJPSI34GeV");
+ TNtuple *nt=(TNtuple*)f->Get("fNtupleJPSI46GeV"); // N-Tuples in 4 - 6 Gev
+ nt->ls();
+
+ aBtoJPSItoEle->ReadCandidates(nt,x,m,n); // read N-Tuples
+ printf("+++\n+++
+ Number of total candidates (prim J/psi + secondary J/psi + bkg) ---> %d candidates
+ \n+++\n",n);
+
+ aBtoJPSItoEle->SetFitHandler(x,m,n); // Set the fit handler with given values of x, m, # of candidates
+
+ aBtoJPSItoEle->CloneMCtemplate(hCsiMCPithya); // clone MC template and copy internally
+ // in this way any model can be setted from outside
+ //aBtoJPSItoEle->CloneMCtemplate(hCsiMCEvtGen); // clone MC template and copy internally
+ // in this way any model can be setted from outside
+
+ aBtoJPSItoEle->SetCsiMC(); // Pass the MC template to the CDF fit function
+
+ AliBtoJPSItoEleCDFfitHandler* fithandler = aBtoJPSItoEle->GetCDFFitHandler(); // Get the fit handler
+
+ //
+ // Set some fit options through the handler class
+ //
+ fithandler->SetErrorDef(0.5); // tells Minuit that the error interval is the one in which
+ // the function differs from the minimum for less than setted value
+
+ fithandler->SetPrintStatus(kTRUE);
+
+ fithandler->SetParamStartValues(paramInputValues);
+ fithandler->SetResolutionConstants();
+ fithandler->SetCrystalBallFunction(kTRUE);
+ fithandler->SetMassWndLow(0.6);
+ fithandler->SetMassWndHigh(0.4);
+ //fithandler->FixParam(5,kTRUE); // Fix some parameters to their input values
+ //fithandler->FixParam(14,kTRUE);
+ //fithandler->FixParam(15,kTRUE);
+ fithandler->SetMassWndLow(0.341696); // ----> M_low = 2.75 GeV
+ fithandler->SetMassWndHigh(0.308304); // ----> M_high = 3.4 GeV
+
+ aBtoJPSItoEle->DoMinimization();
+
+ f->Close();
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff