From f9fd141209ac8c2ae35903007688da239c91f1fd Mon Sep 17 00:00:00 2001 From: dainese Date: Mon, 12 Oct 2009 12:39:01 +0000 Subject: [PATCH] Update of B->Jpsi->ee analysis classes --- PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx | 45 +- PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h | 16 +- .../AliAnalysisTaskSEBtoJPSItoEle.cxx | 151 +++-- .../AliAnalysisTaskSEBtoJPSItoEle.h | 4 +- PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx | 536 +++++++++++------- PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h | 108 ++-- .../AliBtoJPSItoEleCDFfitHandler.cxx | 159 ++++-- .../AliBtoJPSItoEleCDFfitHandler.h | 25 +- PWG3/vertexingHF/macros/FitCDFLocal.C | 110 ++++ 9 files changed, 768 insertions(+), 386 deletions(-) create mode 100644 PWG3/vertexingHF/macros/FitCDFLocal.C diff --git a/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx b/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx index a2ff11868ab..4814efa964f 100644 --- a/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx +++ b/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx @@ -64,41 +64,16 @@ AliAnalysisBtoJPSItoEle::~AliAnalysisBtoJPSItoEle() 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) { // @@ -122,13 +97,27 @@ void AliAnalysisBtoJPSItoEle::ReadCandidates(TNtuple* nt, Double_t* &pseudoprope 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; + + } +} diff --git a/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h b/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h index 9000d4011c5..f432564c84b 100644 --- a/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h +++ b/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h @@ -11,9 +11,11 @@ // Contact: Carmelo.Digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it //------------------------------------------------------------------------- class TNtuple ; -class AliBtoJPSItoEleCDFfitFCN ; -class AliBtoJPSItoEleCDFfitHandler ; +//class AliBtoJPSItoEleCDFfitFCN ; +//class AliBtoJPSItoEleCDFfitHandler ; #include "TH1F.h" +#include "AliBtoJPSItoEleCDFfitHandler.h" +#include "AliBtoJPSItoEleCDFfitFCN.h" //_________________________________________________________________________ class AliAnalysisBtoJPSItoEle : public TNamed { @@ -24,14 +26,16 @@ class AliAnalysisBtoJPSItoEle : public TNamed { AliAnalysisBtoJPSItoEle& operator=(const AliAnalysisBtoJPSItoEle& source); virtual ~AliAnalysisBtoJPSItoEle(); - Int_t DoMinimization(Double_t* x,Double_t* m, Int_t n); + Int_t DoMinimization(); void ReadCandidates(TNtuple* nt, Double_t* &x, Double_t* &m, Int_t& n); // primary JPSI + secondary JPSI + bkg couples void SetPtBin(Int_t BinNum) { fPtBin = BinNum ; } - void SetCsiMC(TH1F* MCtemplate); - void SetResolutionConstants(Int_t BinNum); + void SetCsiMC(); + //void SetResolutionConstants(Int_t BinNum); + void SetFitHandler(Double_t* /*pseudoproper*/, Double_t* /*inv mass*/, Int_t /*candidates*/); void CloneMCtemplate(const TH1F* MCtemplate) {fMCtemplate = (TH1F*)MCtemplate->Clone("fMCtemplate");} Double_t* GetResolutionConstants(Double_t* resolutionConst); + AliBtoJPSItoEleCDFfitHandler* GetCDFFitHandler() const { return fFCNfunction ; } Int_t GetPtBin() const { return fPtBin ; } private: @@ -40,7 +44,7 @@ class AliAnalysisBtoJPSItoEle : public TNamed { Int_t fPtBin; // number of pt bin in which the analysis is performes TH1F* fMCtemplate; //! template of the MC distribution for the x distribution of the secondary J/psi - ClassDef(AliAnalysisBtoJPSItoEle,0); // AliAnalysisBtoJPSItoEle class + ClassDef(AliAnalysisBtoJPSItoEle,1); // AliAnalysisBtoJPSItoEle class }; #endif diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx index 1f35dd6dbfb..66723ed9b7f 100644 --- a/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx +++ b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx @@ -95,11 +95,6 @@ AliAnalysisTaskSEBtoJPSItoEle::~AliAnalysisTaskSEBtoJPSItoEle() 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"); @@ -124,7 +119,7 @@ void AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects() 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); @@ -182,15 +177,24 @@ void AliAnalysisTaskSEBtoJPSItoEle::UserExec(Option_t */*option*/) return; } + // load MC particles + TClonesArray* mcArray = + dynamic_cast(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.}; @@ -202,8 +206,14 @@ void AliAnalysisTaskSEBtoJPSItoEle::UserExec(Option_t */*option*/) } //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()); @@ -225,6 +235,8 @@ void AliAnalysisTaskSEBtoJPSItoEle::UserExec(Option_t */*option*/) fNtupleJPSI->Fill(pseudoProperDecayTime,d->InvMassJPSIee()); // fill N-tuple with X,M values + }// + } // end of JPSItoEle candidates selection according to cuts if(unsetvtx) d->UnsetOwnPrimaryVtx(); @@ -278,7 +290,8 @@ void AliAnalysisTaskSEBtoJPSItoEle::Terminate(Option_t */*option*/) } printf("+++\n+++ MC template histo copied ---> OK\n+++\n"); - fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand); + //fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand); + fCdfFit->DoMinimization(); } return; @@ -300,19 +313,17 @@ void AliAnalysisTaskSEBtoJPSItoEle::SetCutsD0(const Double_t cuts[9]) 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(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 = @@ -321,7 +332,9 @@ void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(const AliAODEvent* aodEvent, c 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 @@ -335,7 +348,19 @@ void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(const AliAODEvent* aodEvent, c 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); @@ -383,36 +408,61 @@ void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(const AliAODEvent* aodEvent, c } } + 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(); @@ -421,11 +471,10 @@ void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(const AliAODEvent* aodEvent, c // 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 } } diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h index c10eb3d59f8..2ebb4f45708 100644 --- a/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h +++ b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h @@ -37,7 +37,7 @@ class AliAnalysisTaskSEBtoJPSItoEle : public AliAnalysisTaskSE void SetPtCuts(const Double_t ptcuts[2]); void SetAODMCInfo(Bool_t OkMCInfo) { fOkAODMC = OkMCInfo;} void SetMinimize(Bool_t OkMinimize) { fOkMinimize = OkMinimize;} - void ReadAODMCInfo(const AliAODEvent* aodEv, const TClonesArray* inArray); + void ReadAODMCInfo(AliAODEvent* aodEv, const TClonesArray* inArray); void SetPathMCfile (TString mcfilename="CsiMCfromKine_10PtBins.root") {fNameMCfile = mcfilename;} TH1F* OpenLocalMCFile(TString dataDir, Int_t nbin); @@ -64,6 +64,6 @@ class AliAnalysisTaskSEBtoJPSItoEle : public AliAnalysisTaskSE Double_t fCuts[9]; // cuts for N-tuple values selection Double_t fPtCuts[2]; // Max and min pt of the candidates - ClassDef(AliAnalysisTaskSEBtoJPSItoEle,0); // AliAnalysisTaskSE for the reconstruction of heavy-flavour decay candidates + ClassDef(AliAnalysisTaskSEBtoJPSItoEle,1); // AliAnalysisTaskSE for the reconstruction of heavy-flavour decay candidates }; #endif diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx index 0f693a06a9d..8dbe41e33d7 100644 --- a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx +++ b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx @@ -32,7 +32,14 @@ AliBtoJPSItoEleCDFfitFCN::AliBtoJPSItoEleCDFfitFCN() : 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.), @@ -41,12 +48,12 @@ fCrystalBallParam(kFALSE) // // 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"); } //_________________________________________________________________________________________________ @@ -56,6 +63,13 @@ fFPlus(source.fFPlus), 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), @@ -64,8 +78,8 @@ fCrystalBallParam(source.fCrystalBallParam) // // 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) @@ -78,11 +92,18 @@ AliBtoJPSItoEleCDFfitFCN& AliBtoJPSItoEleCDFfitFCN::operator=(const AliBtoJPSIto 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; } @@ -90,28 +111,27 @@ AliBtoJPSItoEleCDFfitFCN& AliBtoJPSItoEleCDFfitFCN::operator=(const AliBtoJPSIto 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); @@ -124,45 +144,143 @@ void AliBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parameters) // // 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()); @@ -181,66 +299,43 @@ void AliBtoJPSItoEleCDFfitFCN::PrintStatus() 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 @@ -255,17 +350,19 @@ Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m) c //_________________________________________________________________________________________________ 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; } //_________________________________________________________________________________________________ @@ -275,26 +372,24 @@ Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t m) const // 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]; @@ -306,13 +401,14 @@ Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t m) const //_________________________________________________________________________________________________ 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; @@ -320,29 +416,39 @@ Double_t AliBtoJPSItoEleCDFfitFCN::FunB(Double_t x) const 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; @@ -350,41 +456,50 @@ Double_t AliBtoJPSItoEleCDFfitFCN::CsiMC(Double_t x) const //_________________________________________________________________________________________________ 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; @@ -392,25 +507,30 @@ Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x) const 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; @@ -418,26 +538,29 @@ Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x) const 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; @@ -446,39 +569,50 @@ Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x) const 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 + } + diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h index 297b61d4521..47daa98860d 100644 --- a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h +++ b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h @@ -17,6 +17,7 @@ #include #include "TH1F.h" #include "TMath.h" +#include "TRandom3.h" enum IntegralType {kBkg, kBkgNorm, @@ -39,24 +40,34 @@ class AliBtoJPSItoEleCDFfitFCN : public TNamed { Double_t EvaluateLikelihood(const Double_t* pseudoproperdecaytime, const Double_t* invariantmass, const Int_t ncand); - Double_t GetFPlus() const {return fFPlus;} - Double_t GetFMinus() const {return fFMinus;} - Double_t GetFSym() const {return fFSym;} - Double_t GetRadius() const {return fParameters[0];} - Double_t GetTheta() const {return fParameters[1];} - Double_t GetPhi() const {return fParameters[2];} - Double_t GetLamPlus() const {return fParameters[3];} - Double_t GetLamMinus() const {return fParameters[4];} - Double_t GetLamSym() const {return fParameters[5];} - Double_t GetMassSlope() const {return fParameters[6];} - Double_t GetFractionJpsiFromBeauty() const {return fParameters[7];} - Double_t GetFsig() const {return fParameters[8];} - Double_t GetCrystalBallMmean() const {return fParameters[9];} - Double_t GetCrystalBallNexp() const {return fParameters[10];} - Double_t GetCrystalBallSigma() const {return fParameters[11];} - Double_t GetCrystalBallAlpha() const {return fParameters[12];} - Double_t GetIntegral() const {return fIntegral;} - Bool_t GetCrystalBallParam() const {return fCrystalBallParam;} + Double_t GetFPlus() const { return fFPlus; } + Double_t GetFMinus() const { return fFMinus; } + Double_t GetFSym() const { return fFSym; } + Double_t GetRadius() const { return fParameters[0]; } + Double_t GetTheta() const { return fParameters[1]; } + Double_t GetPhi() const { return fParameters[2]; } + Double_t GetLamPlus() const { return fParameters[3]; } + Double_t GetLamMinus() const { return fParameters[4]; } + Double_t GetLamSym() const { return fParameters[5]; } + Double_t GetMassSlope() const { return fParameters[6]; } + Double_t GetFractionJpsiFromBeauty() const { return fParameters[7]; } + Double_t GetFsig() const { return fParameters[8]; } + Double_t GetCrystalBallMmean() const { return fParameters[9]; } + Double_t GetCrystalBallNexp() const { return fParameters[10]; } + Double_t GetCrystalBallSigma() const { return fParameters[11]; } + Double_t GetCrystalBallAlpha() const { return fParameters[12]; } + Double_t GetCrystalBallNorm() const { return fParameters[13]; } + Double_t GetSigmaResol() const { return fParameters[14]; } + Double_t GetNResol() const { return fParameters[15]; } + Double_t GetIntegral() const { return fIntegral; } + Double_t GetIntegralFunB() const { return fintxFunB; } + Double_t GetIntegralBkgPos() const { return fintxDecayTimeBkgPos; } + Double_t GetIntegralBkgNeg() const { return fintxDecayTimeBkgNeg; } + Double_t GetIntegralBkgSym() const { return fintxDecayTimeBkgSym; } + Double_t GetIntegralMassSig() const { return fintmMassSig; } + Double_t GetIntegralRes() const { return fintxRes; } + Double_t GetIntegralMassBkg() const { return fintmMassBkg; } + Bool_t GetCrystalBallParam() const { return fCrystalBallParam; } void SetFPlus(Double_t plus) {fFPlus = plus;} void SetFMinus(Double_t minus) {fFMinus = minus;} @@ -74,14 +85,26 @@ class AliBtoJPSItoEleCDFfitFCN : public TNamed { void SetCrystalBallNexp(Double_t CrystalBallNexp) {fParameters[10] = CrystalBallNexp;} void SetCrystalBallSigma(Double_t CrystalBallSigma) {fParameters[11] = CrystalBallSigma;} void SetCrystalBallAlpha(Double_t CrystalBallAlpha) {fParameters[12] = CrystalBallAlpha;} + void SetCrystalBallNorm(Double_t CrystalBallNorm) {fParameters[13] = CrystalBallNorm;} + void SetSigmaResol(Double_t SigmaResol) {fParameters[14] = SigmaResol;} + void SetNResol(Double_t NResol) {fParameters[15] = NResol;} void SetAllParameters(const Double_t* parameters); - void SetIntegral(Double_t integral) {fIntegral = integral;} + void SetIntegral(Double_t integral) { fIntegral = integral; } + void SetIntegralFunB(Double_t integral) { fintxFunB = integral; } + void SetIntegralBkgPos(Double_t integral) { fintxDecayTimeBkgPos = integral; } + void SetIntegralBkgNeg(Double_t integral) { fintxDecayTimeBkgNeg = integral; } + void SetIntegralBkgSym(Double_t integral) { fintxDecayTimeBkgSym = integral; } + void SetIntegralMassSig(Double_t integral) { fintmMassSig = integral; } + void SetIntegralRes(Double_t integral) { fintxRes = integral; } + void SetIntegralMassBkg(Double_t integral) { fintmMassBkg = integral; } + void SetCsiMC(const TH1F* MCtemplate) {fhCsiMC = (TH1F*)MCtemplate->Clone("fhCsiMC");} - void SetResolutionConstants(Int_t BinNum); - void SetMassWndHigh(Double_t limit) { fMassWndHigh = TDatabasePDG::Instance()->GetParticle(443)->Mass() + limit ;}//here use pdg code instead - void SetMassWndLow(Double_t limit) { fMassWndLow = TDatabasePDG::Instance()->GetParticle(443)->Mass() - limit ;}//here use pdg code instead - void SetCrystalBallParam(Bool_t okCB) {fCrystalBallParam = okCB;} + + void SetResolutionConstants(); + void SetMassWndHigh(Double_t limit) { fMassWndHigh = TDatabasePDG::Instance()->GetParticle(443)->Mass() + limit ;} + void SetMassWndLow(Double_t limit) { fMassWndLow = TDatabasePDG::Instance()->GetParticle(443)->Mass() - limit ;} + void SetCrystalBallFunction(Bool_t okCB) {fCrystalBallParam = okCB;} void ConvertFromSpherical() { fFPlus = TMath::Power((fParameters[0]*TMath::Cos(fParameters[1])),2.); fFMinus = TMath::Power((fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2])),2.); @@ -95,7 +118,7 @@ class AliBtoJPSItoEleCDFfitFCN : public TNamed { private: // - Double_t fParameters[13]; /* par[0] = fRadius; + Double_t fParameters[16]; /* par[0] = fRadius; par[1] = fTheta; par[2] = fPhi; par[3] = fOneOvLamPlus; @@ -107,18 +130,29 @@ class AliBtoJPSItoEleCDFfitFCN : public TNamed { par[9] = fCrystalBallMmean; par[10] = fCrystalBallNexp; par[11] = fCrystalBallSigma; - par[12] = fCrystalBallAlpha;*/ - - Double_t fFPlus; // parameters of the log-likelihood function - Double_t fFMinus; // Slopes of the x distributions of the background - Double_t fFSym; // functions - - Double_t fIntegral; // integral values of log-likelihood function terms - TH1F *fhCsiMC; // X distribution used as MC template for JPSI from B - Double_t fResolutionConstants[7]; // constants for the parametrized resolution function R(X) - Double_t fMassWndHigh; // JPSI Mass window higher limit - Double_t fMassWndLow; // JPSI Mass window lower limit - Bool_t fCrystalBallParam; // Boolean to switch to Crystall Ball parameterisation + par[12] = fCrystalBallAlpha; + par[13] = fCrystalBallNorm; + par[14] = fSigmaResol; + par[15] = fNResol; */ + + + Double_t fFPlus; // parameters of the log-likelihood function + Double_t fFMinus; // Slopes of the x distributions of the background + Double_t fFSym; // functions + + Double_t fIntegral; // integral values of log-likelihood function terms + Double_t fintxFunB; + Double_t fintxDecayTimeBkgPos; + Double_t fintxDecayTimeBkgNeg; + Double_t fintxDecayTimeBkgSym; + Double_t fintmMassSig; + Double_t fintxRes; + Double_t fintmMassBkg; + TH1F *fhCsiMC; // X distribution used as MC template for JPSI from B + Double_t fResolutionConstants[5]; // constants for the parametrized resolution function R(X) + Double_t fMassWndHigh; // JPSI Mass window higher limit + Double_t fMassWndLow; // JPSI Mass window lower limit + Bool_t fCrystalBallParam; // Boolean to switch to Crystall Ball parameterisation //// @@ -156,7 +190,7 @@ class AliBtoJPSItoEleCDFfitFCN : public TNamed { Double_t ResolutionFunc(Double_t x) const ; // - ClassDef (AliBtoJPSItoEleCDFfitFCN,0); // Unbinned log-likelihood fit + ClassDef (AliBtoJPSItoEleCDFfitFCN,1); // Unbinned log-likelihood fit }; diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx index 431e1509152..6c7851444c5 100644 --- a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx +++ b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx @@ -43,6 +43,8 @@ void CDFFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t i //_________________________________________________________________________________________________ AliBtoJPSItoEleCDFfitHandler::AliBtoJPSItoEleCDFfitHandler(): +fIsParamFixed(16), +fPrintStatus(kFALSE), fUp(0), fX(0x0), fM(0x0), @@ -56,6 +58,8 @@ fNcand(0) //_________________________________________________________________________________________________ AliBtoJPSItoEleCDFfitHandler::AliBtoJPSItoEleCDFfitHandler(Double_t* decaytime, Double_t* invariantmass, Int_t ncand) : +fIsParamFixed(16), +fPrintStatus(kFALSE), fUp(0), fX(decaytime), fM(invariantmass), @@ -66,10 +70,24 @@ fNcand(ncand) // 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)); } //___________________________________________________________________________ @@ -79,18 +97,22 @@ AliBtoJPSItoEleCDFfitHandler& AliBtoJPSItoEleCDFfitHandler::operator=(const AliB // 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), @@ -101,7 +123,7 @@ fNcand(c.fNcand) // Copy Constructor // } -//________________________________________________________________________ +//_______________________________________________________________________________________ AliBtoJPSItoEleCDFfitHandler::~AliBtoJPSItoEleCDFfitHandler() { // @@ -109,58 +131,40 @@ 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); @@ -178,8 +182,7 @@ void AliBtoJPSItoEleCDFfitHandler::CdfFCN(Int_t & /* npar */, fLikely->SetAllParameters(par); fLikely->ConvertFromSpherical(); fLikely->ComputeIntegral(); - //printf("\n+++\n+++\n+++\n"); - //fLikely->PrintStatus(); + if(fPrintStatus)fLikely->PrintStatus(); TStopwatch t; t.Start(); @@ -193,3 +196,43 @@ void AliBtoJPSItoEleCDFfitHandler::CdfFCN(Int_t & /* npar */, 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); +} + diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h index 534934be30f..9c55f255e4a 100644 --- a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h +++ b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h @@ -12,9 +12,12 @@ //____________________________________________________________________________ #include +#include #include "TMath.h" #include "AliBtoJPSItoEleCDFfitFCN.h" -//#include "AliLog.h" +#include "AliLog.h" + +class TBits; class AliBtoJPSItoEleCDFfitHandler : public TNamed { public: @@ -24,8 +27,21 @@ class AliBtoJPSItoEleCDFfitHandler : public TNamed { AliBtoJPSItoEleCDFfitHandler(const AliBtoJPSItoEleCDFfitHandler& c); AliBtoJPSItoEleCDFfitHandler(Double_t* decaytime, Double_t* invariantmass, Int_t ncand); ~AliBtoJPSItoEleCDFfitHandler(); - Double_t Up() const { return fUp*TMath::Sqrt(fLikely->GetIntegral()); } + //Double_t Up() const { return fUp*TMath::Sqrt(fLikely->GetIntegral()); } + Double_t Up() const { return fUp; } void SetErrorDef(Double_t up) {fUp = up;} + void SetPrintStatus(Bool_t okPrintStatus) { fPrintStatus = okPrintStatus; } + Bool_t GetPrintStatus() { return fPrintStatus ; } + void SetParamStartValues (Double_t*); + Double_t* GetStartParamValues() { return fParamStartValues; } + TBits GetFixedParamList() { return fIsParamFixed; } + void FixParam(UInt_t param, Bool_t value) { fIsParamFixed.SetBitNumber(param,value); } + void FixAllParam(Bool_t value) { for(UInt_t par=0;par<16;par++) fIsParamFixed.SetBitNumber(par,value); } + Bool_t IsParamFixed(UInt_t param) { return fIsParamFixed.TestBitNumber(param); } + void SetResolutionConstants(); + void SetCrystalBallFunction(Bool_t okCB); + void SetMassWndHigh(Double_t limit); + void SetMassWndLow(Double_t limit); Double_t operator()(const Double_t* par) const ; void CdfFCN(Int_t & /* npar */, Double_t * /* gin */, Double_t &f, Double_t *par, Int_t /* iflag */); @@ -37,13 +53,16 @@ class AliBtoJPSItoEleCDFfitHandler : public TNamed { private: // + TBits fIsParamFixed; //array of bits: 0 = param free; 1 = param fixed; + Bool_t fPrintStatus; //flag to enable the prit out of the algorithm at each step + Double_t fParamStartValues[16]; //array of parameters input value Double_t fUp; //error definition Double_t* fX; //pseudo-proper decay time X Double_t* fM; //invariant mass M AliBtoJPSItoEleCDFfitFCN* fLikely; //Log likelihood function Int_t fNcand; //number of candidates // - ClassDef(AliBtoJPSItoEleCDFfitHandler,0); + ClassDef(AliBtoJPSItoEleCDFfitHandler,1); }; #endif diff --git a/PWG3/vertexingHF/macros/FitCDFLocal.C b/PWG3/vertexingHF/macros/FitCDFLocal.C new file mode 100644 index 00000000000..851e14a9e9e --- /dev/null +++ b/PWG3/vertexingHF/macros/FitCDFLocal.C @@ -0,0 +1,110 @@ +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(); +} -- 2.43.0