From 7b54a4836407cad3f023c261cdc2e3327bcefe68 Mon Sep 17 00:00:00 2001 From: dainese Date: Thu, 16 Apr 2009 09:17:30 +0000 Subject: [PATCH] New set of classes for B->J/psi->ee analysis and fit a la CDF (Carmelo, Giuseppe) --- PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx | 134 +++++ PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h | 46 ++ .../AliAnalysisTaskSEBtoJPSItoEle.cxx | 459 +++++++++++++++++ .../AliAnalysisTaskSEBtoJPSItoEle.h | 69 +++ .../AliAnalysisTaskSEBtoJPSItoEleTest.C | 94 ++++ PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx | 484 ++++++++++++++++++ PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h | 163 ++++++ .../AliBtoJPSItoEleCDFfitHandler.cxx | 195 +++++++ .../AliBtoJPSItoEleCDFfitHandler.h | 49 ++ 9 files changed, 1693 insertions(+) create mode 100644 PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx create mode 100644 PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h create mode 100644 PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx create mode 100644 PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h create mode 100644 PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEleTest.C create mode 100644 PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx create mode 100644 PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h create mode 100644 PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx create mode 100644 PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h diff --git a/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx b/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx new file mode 100644 index 00000000000..a2ff11868ab --- /dev/null +++ b/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx @@ -0,0 +1,134 @@ +/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//------------------------------------------------------------------------- +// Class AliAnalysisBtoJPSItoEle +// Unbinned log-likelihood fit analysis class +// +// Origin: C.Di Giglio +// Contact: Carmelo.Digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it +//------------------------------------------------------------------------- +class TH1F ; +#include "TNtuple.h" +#include "TMath.h" + +#include "AliBtoJPSItoEleCDFfitFCN.h" +#include "AliBtoJPSItoEleCDFfitHandler.h" +#include "AliAnalysisBtoJPSItoEle.h" +#include "AliLog.h" + +ClassImp(AliAnalysisBtoJPSItoEle) + +//_______________________________________________________________________________ +AliAnalysisBtoJPSItoEle::AliAnalysisBtoJPSItoEle() : +fFCNfunction(0), +fPtBin(0), +fMCtemplate(0) +{ + // + // default constructor + // +} +//___________________________________________________________________________________ +AliAnalysisBtoJPSItoEle::AliAnalysisBtoJPSItoEle(const AliAnalysisBtoJPSItoEle& source) : +TNamed(source), +fFCNfunction(source.fFCNfunction), +fPtBin(source.fPtBin), +fMCtemplate(source.fMCtemplate) +{ + // + // copy constructor + // +} +//_________________________________________________________________________________________________ + +AliAnalysisBtoJPSItoEle &AliAnalysisBtoJPSItoEle::operator=(const AliAnalysisBtoJPSItoEle& source) +{ + // + // assignment operator + // + if(&source == this) return *this; + fFCNfunction = source.fFCNfunction; + fPtBin = source.fPtBin; + fMCtemplate = source.fMCtemplate; + + return *this; +} +//_________________________________________________________________________________________________ +AliAnalysisBtoJPSItoEle::~AliAnalysisBtoJPSItoEle() +{ + // + // destructor + // + delete fFCNfunction; + delete fMCtemplate; +} +//_________________________________________________________________________________________________ +Int_t AliAnalysisBtoJPSItoEle::DoMinimization(Double_t* x, + Double_t* m, Int_t ncand) +{ + // + // 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) +{ + // + // Read N-tuple with X and M values + // + Float_t mJPSI = 0; Float_t x = 0; + Int_t nentries = 0; + ncand=0; + nt->SetBranchAddress("Mass",&mJPSI); + nt->SetBranchAddress("Xdecaytime",&x); + nentries = (Int_t)nt->GetEntries(); + pseudoproper = new Double_t[nentries]; + invmass = new Double_t[nentries]; + for(Int_t i = 0; i < nentries; i++) { + nt->GetEntry(i); + ncand++; + pseudoproper[i]=(Double_t)(10000*x); + invmass[i]=(Double_t)mJPSI; + } + + return; +} +//_________________________________________________________________________________________________ +void AliAnalysisBtoJPSItoEle::SetCsiMC(TH1F* MCtemplates) +{ + // + // Sets X distribution used as MC template for JPSI from B + // + fFCNfunction->LikelihoodPointer()->SetCsiMC(MCtemplates); + + return; +} +//_________________________________________________________________________________________________ diff --git a/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h b/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h new file mode 100644 index 00000000000..9000d4011c5 --- /dev/null +++ b/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h @@ -0,0 +1,46 @@ +#ifndef ALIANALYSISBTOJPSITOELE_H +#define ALIANALYSISBTOJPSITOELE_H +/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//------------------------------------------------------------------------- +// Class AliAnalysisBtoJPSItoEle +// Unbinned log-likelihood fit analysis class +// +// Origin: C.Di Giglio +// Contact: Carmelo.Digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it +//------------------------------------------------------------------------- +class TNtuple ; +class AliBtoJPSItoEleCDFfitFCN ; +class AliBtoJPSItoEleCDFfitHandler ; +#include "TH1F.h" + +//_________________________________________________________________________ +class AliAnalysisBtoJPSItoEle : public TNamed { + public: + // + AliAnalysisBtoJPSItoEle(); + AliAnalysisBtoJPSItoEle(const AliAnalysisBtoJPSItoEle& source); + AliAnalysisBtoJPSItoEle& operator=(const AliAnalysisBtoJPSItoEle& source); + virtual ~AliAnalysisBtoJPSItoEle(); + + Int_t DoMinimization(Double_t* x,Double_t* m, Int_t n); + 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 CloneMCtemplate(const TH1F* MCtemplate) {fMCtemplate = (TH1F*)MCtemplate->Clone("fMCtemplate");} + Double_t* GetResolutionConstants(Double_t* resolutionConst); + Int_t GetPtBin() const { return fPtBin ; } + + private: + // + AliBtoJPSItoEleCDFfitHandler* fFCNfunction; //! pointer to the interface class + 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 +}; + +#endif diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx new file mode 100644 index 00000000000..70310fb6e91 --- /dev/null +++ b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx @@ -0,0 +1,459 @@ +/************************************************************************** + * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ +//************************************************************************* +// Class AliAnalysisTaskSEBtoJPSItoEle +// AliAnalysisTaskSE for the reconstruction +// of heavy-flavour decay candidates +// Author: C.Di Giglio, carmelo.digiglio@ba.infn.it +//************************************************************************* +class AliAnalysisTaskSE; +class AliESDEvent; +class AliVEvent; +class iostream; +class TSystem; +class TROOT; +#include +#include +#include + +#include "AliAODEvent.h" +#include "AliAODMCParticle.h" +#include "AliAODMCHeader.h" +#include "AliAODRecoDecayHF2Prong.h" + +#include "AliAnalysisBtoJPSItoEle.h" +#include "AliAnalysisTaskSEBtoJPSItoEle.h" + +ClassImp(AliAnalysisTaskSEBtoJPSItoEle) + +//______________________________________________________________________________ +AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle() : +AliAnalysisTaskSE(), +fOutput(0), +fCdfFit(0), +fNtupleJPSI(0), +fhDecayTimeMCjpsifromB(0), +fhDecayTime(0), +fhInvMass(0), +fhD0(0), +fhD0D0(0), +fhCosThetaStar(0), +fhCosThetaPointing(0), +fhCtsVsD0D0(0), +fOkAODMC(kFALSE), +fNameMCfile(""), +fOkMinimize(kFALSE) +{ + // default constructor + + // Output slot #1 writes into a TList container + DefineOutput(1,TList::Class()); //My private output +} +//_________________________________________________________________________________________________ +AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle(const char* name) : +AliAnalysisTaskSE(name), +fOutput(0), +fCdfFit(0), +fNtupleJPSI(0), +fhDecayTimeMCjpsifromB(0), +fhDecayTime(0), +fhInvMass(0), +fhD0(0), +fhD0D0(0), +fhCosThetaStar(0), +fhCosThetaPointing(0), +fhCtsVsD0D0(0), +fOkAODMC(kFALSE), +fNameMCfile(""), +fOkMinimize(kFALSE) +{ + // default constructor + + // Output slot #1 writes into a TList container + DefineOutput(1,TList::Class()); //My private output +} +//_________________________________________________________________________________________________ +AliAnalysisTaskSEBtoJPSItoEle::~AliAnalysisTaskSEBtoJPSItoEle() +{ + // destructor + + if (fOutput) { + delete fOutput; + fOutput = 0; + } +} +//_________________________________________________________________________________________________ +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"); + + return; + +} +//_________________________________________________________________________________________________ +void AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects() +{ + // Create the output container + + if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects() \n"); + + // Several histograms are more conveniently managed in a TList + fOutput = new TList(); + fOutput->SetOwner(); + + if(fOkAODMC){ + fhDecayTimeMCjpsifromB = new TH1F("fhDecayTimeMCjpsifromB", "Secondary J/#Psi Monte Carlo pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.); + fhDecayTimeMCjpsifromB->Sumw2(); + fhDecayTimeMCjpsifromB->SetMinimum(0); + fOutput->Add(fhDecayTimeMCjpsifromB); + } + + fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.); + fhDecayTime->Sumw2(); + fhDecayTime->SetMinimum(0); + fOutput->Add(fhDecayTime); + + fhInvMass = new TH1F("fhInvMass", "J/#Psi invariant mass; M [GeV]; Entries",100,0.,3.2); + fhInvMass->Sumw2(); + fhInvMass->SetMinimum(0); + fOutput->Add(fhInvMass); + + fhD0 = new TH1F("fhD0", "Impact parameter; D0 [#mu m]; Entries",100,-5000.,5000.); + fhD0->Sumw2(); + fhD0->SetMinimum(0); + fOutput->Add(fhD0); + + fhD0D0 = new TH1F("fhD0D0", "Product of impact parameters; D0D0 [#mu m^2]; Entries",100,-100000.,100000.); + fhD0D0->Sumw2(); + fhD0D0->SetMinimum(0); + fOutput->Add(fhD0D0); + + fhCosThetaStar = new TH1F("fhCosThetaStar", "Cosine of decay angle; Cosine Theta Star; Entries",50,-1.2,1.2); + fhCosThetaStar->Sumw2(); + fhCosThetaStar->SetMinimum(0); + fOutput->Add(fhCosThetaStar); + + fhCosThetaPointing = new TH1F("fhCosThetaPointing", "Cosine of pointing angle; Cosine Pointing Angle; Entries",100,-1,1); + fhCosThetaPointing->Sumw2(); + fhCosThetaPointing->SetMinimum(0); + fOutput->Add(fhCosThetaPointing); + + fhCtsVsD0D0 = new TH2F("fhCtsVsD0D0", "Cosine theta star Vs. D0; Cosine theta star; D0D0",50,-1,1,100,-100000,100000); + fhCtsVsD0D0->Sumw2(); + fhCtsVsD0D0->SetMinimum(0); + fOutput->Add(fhCtsVsD0D0); + + fNtupleJPSI = new TNtuple("fNtupleJPSI","J/#Psi pseudo-proper decay time & invariant mass","Xdecaytime:Mass"); + fOutput->Add(fNtupleJPSI); + + return; +} +//_________________________________________________________________________________________________ +void AliAnalysisTaskSEBtoJPSItoEle::UserExec(Option_t */*option*/) +{ + // Execute analysis for current event: + + AliAODEvent *aod = dynamic_cast (InputEvent()); + + // AOD primary vertex + AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex(); + + // load JPSI candidates + TClonesArray *inputArrayJPSItoEle = + (TClonesArray*)aod->GetList()->FindObject("JPSItoEle"); + if(!inputArrayJPSItoEle) { + printf("AliAnalysisTaskSECompareHF::UserExec: JPSItoEle branch not found!\n"); + return; + } + + // Read AOD Monte Carlo information + if(fOkAODMC) ReadAODMCInfo(aod,inputArrayJPSItoEle); + + // loop over J/Psi->ee candidates + Int_t nInJPSItoEle = inputArrayJPSItoEle->GetEntriesFast(); + printf("Number of B->JPSI->e+e-: %d\n",nInJPSItoEle); + + for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) { + AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayJPSItoEle->UncheckedAt(iJPSItoEle); + //Secondary vertex + Double_t vtxSec[3] = {0.,0.,0.}; + Double_t vtxPrim[3] = {0.,0.,0.}; + d->GetSecondaryVtx(vtxSec); + Bool_t unsetvtx=kFALSE; + if(!d->GetOwnPrimaryVtx()) { + d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables + unsetvtx=kTRUE; + } + //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 + fhInvMass->Fill(d->InvMassJPSIee()); + fhD0->Fill(10000*d->ImpParXY()); + fhD0D0->Fill(1e8*d->Prodd0d0()); + fhCosThetaStar->Fill(d->CosThetaStarJPSI()); + fhCtsVsD0D0->Fill(d->CosThetaStarJPSI(),1e8*d->Prodd0d0()); + fhCosThetaPointing->Fill(d->CosPointingAngle()); + + // compute X variable + AliAODVertex* primVertex = d->GetOwnPrimaryVtx(); + vtxPrim[0] = primVertex->GetX(); + vtxPrim[1] = primVertex->GetY(); + vtxPrim[2] = primVertex->GetZ(); + Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(d->Px()) + (vtxSec[1]-vtxPrim[1])*(d->Py()))/d->Pt(); + Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/d->Pt(); + + fhDecayTime->Fill(10000*pseudoProperDecayTime); + // Post the data already here + PostData(1,fOutput); + + fNtupleJPSI->Fill(pseudoProperDecayTime,d->InvMassJPSIee()); // fill N-tuple with X,M values + + } // end of JPSItoEle candidates selection according to cuts + + if(unsetvtx) d->UnsetOwnPrimaryVtx(); + + }// end loop on JPSI to ele candidates + +} +//_________________________________________________________________________________________________ +void AliAnalysisTaskSEBtoJPSItoEle::Terminate(Option_t */*option*/) +{ + // + // Terminate analysis + // + + if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle: Terminate() \n"); + + fOutput = dynamic_cast (GetOutputData(1)); + if (!fOutput) { + printf("ERROR: fOutput not available\n"); + return; + } + + if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast(fOutput->FindObject("fhDecayTimeMCjpsifromB")); + fhDecayTime = dynamic_cast(fOutput->FindObject("fhDecayTime")); + fhInvMass = dynamic_cast(fOutput->FindObject("fhInvMass")); + fhD0 = dynamic_cast(fOutput->FindObject("fhD0")); + fhD0D0 = dynamic_cast(fOutput->FindObject("fhD0D0")); + fhCosThetaStar = dynamic_cast(fOutput->FindObject("fhCosThetaStar")); + fhCosThetaPointing = dynamic_cast(fOutput->FindObject("fhCosThetaPointing")); + fhCtsVsD0D0 = dynamic_cast(fOutput->FindObject("fhCtsVsD0D0")); + fNtupleJPSI = dynamic_cast(fOutput->FindObject("fNtupleJPSI")); + + if(fOkMinimize){ + Double_t* pseudoProperValues=0x0; + Double_t* invMassValues=0x0; + Int_t ncand=0; + fCdfFit = new AliAnalysisBtoJPSItoEle(); + + printf("+++\n+++ AliAnalysisBtoJPSItoEle object instantiated ---> OK\n+++\n"); + fCdfFit->ReadCandidates(fNtupleJPSI,pseudoProperValues,invMassValues,ncand); + printf("+++\n+++ X and M vectors filled starting from N-tuple ---> OK\n+++\n"); + printf("+++\n+++ Number of candidates ---> %d J/psi \n+++\n",ncand); + fCdfFit->SetPtBin(0); + printf("+++\n+++ Pt bin setted n. ---> %d \n+++\n",fCdfFit->GetPtBin()); + + if(fOkAODMC) {fCdfFit->CloneMCtemplate(fhDecayTimeMCjpsifromB);} + else { + //SetPathMCfile("."); + TH1F* hLocal = OpenLocalMCFile(fNameMCfile,fCdfFit->GetPtBin()); + fCdfFit->CloneMCtemplate(hLocal); + } + + printf("+++\n+++ MC template histo copied ---> OK\n+++\n"); + fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand); + } + + return; +} +//_________________________________________________________________________________________________ +void AliAnalysisTaskSEBtoJPSItoEle::SetPtCuts(const Double_t ptCuts[2]) +{ + // + // apply Pt cuts (lower cut, upper cut) + // + for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i]; +} +//_________________________________________________________________________________________________ +void AliAnalysisTaskSEBtoJPSItoEle::SetCutsD0(const Double_t cuts[9]) +{ + // + // apply D0 and JPSI cuts + // + for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i]; +} +//_________________________________________________________________________________________________ +void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(const 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; + } + + // load MC header + AliAODMCHeader* mcHeader = + (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()); + if(!mcHeader){ + printf("AliAnalysisTaskSEBtoJPSItoEle::UserExec: MC AOD header branch not found!\n"); + } + + Double_t rmax = 0.000005; + Double_t mcPrimVtx[3]; + + // get MC primary Vtx + mcHeader->GetVertex(mcPrimVtx); + + Int_t nInCandidates = inputArray->GetEntriesFast(); + printf("Number of Candidates for MC analysis: %d\n",nInCandidates); + + Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother; + Int_t labJPSIdaugh0=0; + Int_t labJPSIdaugh1=0; + + for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) { + AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate); + if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only + + AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0); + AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1); + lab0 = trk0->GetLabel(); + lab1 = trk1->GetLabel(); + + AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0); + if(!part0) { + printf("no MC particle\n"); + continue; + } + + while(part0->GetMother()>=0) { + labMother=part0->GetMother(); + part0 = (AliAODMCParticle*)mcArray->At(labMother); + if(!part0) { + printf("no MC mother particle\n"); + break; + } + pdgMother = TMath::Abs(part0->GetPdgCode()); + if(pdgMother==443) {//this for JPSI + labJPSIdaugh0=labMother; + break; + } + } + + AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1); + if(!part1) { + printf("no MC particle\n"); + continue; + } + + while(part1->GetMother()>=0) { + labMother=part1->GetMother(); + part1 = (AliAODMCParticle*)mcArray->At(labMother); + if(!part1) { + printf("no MC mother particle\n"); + break; + } + pdgMother = TMath::Abs(part1->GetPdgCode()); + if(pdgMother==443) {//this for JPSI + labJPSIdaugh1=labMother; + break; + } + } + + 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 + + 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(); + fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1); + + // 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 + } + } + +} +//_________________________________________________________________________________________________ +TH1F* AliAnalysisTaskSEBtoJPSItoEle::OpenLocalMCFile(TString datadir, Int_t binNum) +{ + // + // Open a local file with MC x distribution for JPSI from B-hadron + // + + TH1F* hCsiMC = new TH1F(); + TString useFile = datadir.Data(); + useFile.Append("/CsiMCfromKine_10PtBins.root"); + TFile *f = new TFile(useFile); + Double_t scale = 0.; + char processCsiMCXHisto[1024]; + if(binNum == 0) sprintf(processCsiMCXHisto,"hPseudoProper"); + else sprintf(processCsiMCXHisto,"hPseudoProper%d",binNum); + hCsiMC = (TH1F*)f->Get(processCsiMCXHisto); + scale=1/hCsiMC->Integral(); + hCsiMC->Scale(scale); + printf ("Opening Histo with Csi_MC(x) template with n. %f entries", hCsiMC->GetEntries()); + + return hCsiMC; + +} + diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h new file mode 100644 index 00000000000..f81c8882af0 --- /dev/null +++ b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h @@ -0,0 +1,69 @@ +#ifndef ALIANALYSISTASKSEBTOJPSITOELE_H +#define ALIANALYSISTASKSEBTOJPSITOELE_H +/* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//************************************************************************* +// Class AliAnalysisTaskSEBtoJPSItoEle +// AliAnalysisTaskSE for the reconstruction +// of heavy-flavour decay candidates +// Author: C.Di Giglio, carmelo.digiglio@ba.infn.it +//************************************************************************* + +class TH1F; +class TList; +class AliAnalysisBtoJPSItoEle; +#include +#include +#include +#include "AliAnalysisTaskSE.h" + +class AliAnalysisTaskSEBtoJPSItoEle : public AliAnalysisTaskSE +{ + public: + + AliAnalysisTaskSEBtoJPSItoEle(); + AliAnalysisTaskSEBtoJPSItoEle(const char *name); + virtual ~AliAnalysisTaskSEBtoJPSItoEle(); + + // Implementation of interface methods + virtual void UserCreateOutputObjects(); + virtual void Init(); + virtual void LocalInit() {Init();} + virtual void UserExec(Option_t *option); + virtual void Terminate(Option_t *option); + + void SetCutsD0(const Double_t cutsD0[9]); + 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 SetPathMCfile (TString mcfilename="CsiMCfromKine_10PtBins.root") {fNameMCfile = mcfilename;} + TH1F* OpenLocalMCFile(TString dataDir, Int_t nbin); + + private: + + AliAnalysisTaskSEBtoJPSItoEle(const AliAnalysisTaskSEBtoJPSItoEle &source); + AliAnalysisTaskSEBtoJPSItoEle& operator=(const AliAnalysisTaskSEBtoJPSItoEle& source); + // + TList *fOutput; //! list send on output slot 0 + AliAnalysisBtoJPSItoEle *fCdfFit; // Unbinned log-likelihood minimizer + TNtuple *fNtupleJPSI; // Ntuple of pseudo-proper decay time & invariant mass values + TH1F *fhDecayTimeMCjpsifromB; // Pseudo-proper decay time distribution used as template for JPSIs from B + TH1F *fhDecayTime; // Pseudo-proper decay time distribution + TH1F *fhInvMass; // Invariant mass distribution + TH1F *fhD0; // Impact parameter distribution + TH1F *fhD0D0; // Product of impact parameters distributions + TH1F *fhCosThetaStar; // Cosine of decay angle distribution + TH1F *fhCosThetaPointing; // Cosine of pointing angle distribution + TH2F *fhCtsVsD0D0; // Cos theta star Vs. D0D0 distribution + Bool_t fOkAODMC; // Flag to read AOD monte carlo information + TString fNameMCfile; // Name of the MC file with X template + Bool_t fOkMinimize; // Flag to enable unbinned log-likelihood minimization + + 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 +}; +#endif diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEleTest.C b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEleTest.C new file mode 100644 index 00000000000..60f188b9329 --- /dev/null +++ b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEleTest.C @@ -0,0 +1,94 @@ +void AliAnalysisTaskSEBtoJPSItoEleTest() +{ + // + // Test macro for the AliAnalysisTaskSEBtoJPSItoEle + // starting from AliAOD.root file with HF + Like Sign candidates. + // C.Di Giglio, carmelo.digiglio@ba.infn.it + // + + TString mode="local"; // otherwise, "grid" + //TString mode="grid"; // needs definition of a proper Grid handler + + // load proper libraries + gSystem->Load("libTree.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libPhysics.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libSTEERBase.so"); + gSystem->Load("libESD.so"); + gSystem->Load("libAOD.so"); + gSystem->Load("libANALYSIS.so"); + gSystem->Load("libANALYSISalice.so"); + gSystem->Load("libCORRFW.so"); + gSystem->Load("libPWG3base.so"); + gSystem->Load("libPWG3vertexingHF.so"); + + TChain *chain = 0; + + if(mode=="local") { + // Local files + TString treeName,fileName; + treeName="aodTree"; + fileName="AliAOD.root"; + chain = new TChain(treeName.Data()); + chain->Add(fileName.Data()); + } else if (mode=="grid") { + // Fetch files with AliEn : + const char *collectionfile = "Collection.xml"; + TGrid::Connect("alien://") ; + TAlienCollection *coll = TAlienCollection::Open(collectionfile); + chain = new TChain("aodTree"); + while(coll->Next()) chain->Add(coll->GetTURL("")); + } else { + printf("ERROR: mode has to be \"local\" or \"grid\" \n"); + return; + } + + // Create the analysis manager + AliAnalysisManager *mgr = new AliAnalysisManager("My Manager","My Manager"); + mgr->SetDebugLevel(10); + + // Input + AliAODInputHandler *inputHandler = new AliAODInputHandler(); + inputHandler->AddFriend("AliAOD.VertexingHF.root"); + mgr->SetInputEventHandler(inputHandler); + + // Cdf unbinned log-likelihood fit analysis task + AliAnalysisTaskSEBtoJPSItoEle *hfTask = new AliAnalysisTaskSEBtoJPSItoEle("CdfFitAnalysis"); + hfTask->SetDebugLevel(2); + + mgr->AddTask(hfTask); + + // + // Create containers for input/output + mgr->ConnectInput(hfTask,0,mgr->GetCommonInputContainer()); + + // before v4-17-Release + /*AliAnalysisDataContainer *cinput1 = mgr->CreateContainer("cchain",TChain::Class(), + AliAnalysisManager::kInputContainer); + AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("tree", TTree::Class(), + AliAnalysisManager::kOutputContainer, + "default"); + + mgr->ConnectInput(hfTask,0,cinput1); + mgr->ConnectOutput(hfTask,0,coutput1);*/ + + AliAnalysisDataContainer *coutput = mgr->CreateContainer("coutput",TList::Class(), + AliAnalysisManager::kOutputContainer, + "CdfFit.root"); + + mgr->ConnectOutput(hfTask,1,coutput); + + // + // Run the analysis + // + printf("CHAIN HAS %d ENTRIES\n",(Int_t)chain->GetEntries()); + if(!mgr->InitAnalysis()) return; + + mgr->PrintStatus(); + + mgr->StartAnalysis(mode.Data(),chain); + + return; + +} diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx new file mode 100644 index 00000000000..0f693a06a9d --- /dev/null +++ b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx @@ -0,0 +1,484 @@ +/************************************************************************** + * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ +#include "AliLog.h" +#include "AliBtoJPSItoEleCDFfitFCN.h" + +//_________________________________________________________________________ +// Class AliBtoJPSItoEleCDFfitFCN +// Definition of main function used in +// unbinned log-likelihood fit for +// the channel B -> JPsi + X -> e+e- + X +// +// Origin: C.Di Giglio +// Contact: Carmelo.Digiglio@ba.infn.it , Giuseppe.Bruno@ba.infn.it +//_________________________________________________________________________ + +ClassImp(AliBtoJPSItoEleCDFfitFCN) + +//_________________________________________________________________________________________________ +AliBtoJPSItoEleCDFfitFCN::AliBtoJPSItoEleCDFfitFCN() : +fFPlus(0.), +fFMinus(0.), +fFSym(0.), +fIntegral(0.), +fhCsiMC(0x0), +fMassWndHigh(0.), +fMassWndLow(0.), +fCrystalBallParam(kFALSE) +{ + // + // constructor + // + SetCrystalBallParam(kFALSE); + SetMassWndHigh(0.2); + SetMassWndLow(0.5); + for(Int_t iPar = 0; iPar < 13; iPar++) fParameters[iPar] = 0.; + fParameters[9] = 1.;fParameters[11] = 1.;fParameters[12] = 1.; + for(Int_t index=0; index<7; index++) fResolutionConstants[index] = 0.; + AliInfo("Instance of AliBtoJPSItoEleCDFfitFCN-class created"); +} +//_________________________________________________________________________________________________ +AliBtoJPSItoEleCDFfitFCN::AliBtoJPSItoEleCDFfitFCN(const AliBtoJPSItoEleCDFfitFCN& source) : +TNamed(source), +fFPlus(source.fFPlus), +fFMinus(source.fFMinus), +fFSym(source.fFSym), +fIntegral(source.fIntegral), +fhCsiMC(source.fhCsiMC), +fMassWndHigh(source.fMassWndHigh), +fMassWndLow(source.fMassWndLow), +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]; +} +//_________________________________________________________________________________________________ +AliBtoJPSItoEleCDFfitFCN& AliBtoJPSItoEleCDFfitFCN::operator=(const AliBtoJPSItoEleCDFfitFCN& source) +{ + // + // Assignment operator + // + if(&source == this) return *this; + fFPlus = source.fFPlus; + fFMinus = source.fFMinus; + fFSym = source.fFSym; + fIntegral = source.fIntegral; + 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]; + + return *this; +} +//_________________________________________________________________________________________________ +AliBtoJPSItoEleCDFfitFCN::~AliBtoJPSItoEleCDFfitFCN() +{ + // + // 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.; +} +//_________________________________________________________________________________________________ +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) +// + 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); + } + return ret; +} +//_________________________________________________________________________________________________ +void AliBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parameters) +{ + // + // Sets array of FCN parameters + // + for(Int_t index = 0; index < 13; 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; + 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; + 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); +} +//_________________________________________________________________________________________________ +void AliBtoJPSItoEleCDFfitFCN::PrintStatus() +{ +// +// 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 fPhi ------------------------------------------>> | %f \n", GetPhi()); + printf("actual value of fFPlus ---------------------------------------->> | %f \n", GetFPlus()); + printf("actual value of fFMinus --------------------------------------->> | %f \n", GetFMinus()); + printf("actual value of fFSym ----------------------------------------->> | %f \n", GetFSym()); + printf("actual value of fOneOvLamPlus --------------------------------->> | %f \n", GetLamPlus()); + printf("actual value of fOneOvLamMinus -------------------------------->> | %f \n", GetLamMinus()); + printf("actual value of fOneOvLamSym ---------------------------------->> | %f \n", GetLamSym()); + printf("actual value of fMassBkgSlope --------------------------------->> | %f \n", GetMassSlope()); + printf("actual value of fFractionJpsiFromBeauty ----------------------->> | %f \n", GetFractionJpsiFromBeauty()); + printf("actual value of fFsig ----------------------------------------->> | %f \n", GetFsig()); + if(fCrystalBallParam){ + printf("actual value of fCrystalBallMmean ----------------------------->> | %f \n", GetCrystalBallMmean()); + printf("actual value of fCrystalBallNexp ------------------------------>> | %f \n", GetCrystalBallNexp()); + printf("actual value of fCrystalBallSigma ----------------------------->> | %f \n", GetCrystalBallSigma()); + printf("actual value of fCrystalBallAlpha ----------------------------->> | %f \n", GetCrystalBallAlpha()); + }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("\n"); + printf("Actual value of normalization integral for FCN ---------------->> | %f \n", GetIntegral()); + printf("\n"); +} +//_________________________________________________________________________________________________ +void AliBtoJPSItoEleCDFfitFCN::SetResolutionConstants(Int_t BinNum) +{ +// 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; + } +} +//_________________________________________________________________________________________________ +Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m) const +{ + return fParameters[8]*EvaluateCDFfuncSignalPart(x,m) + (1. - fParameters[8])*EvaluateCDFfuncBkgPart(x,m); +} +//_________________________________________________________________________________________________ +Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m) const +{ + return EvaluateCDFfunc(x,m)/fIntegral; +} +//_________________________________________________________________________________________________ +Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m) const +{ + return EvaluateCDFDecayTimeSigDistr(x)*EvaluateCDFInvMassSigDistr(m); +} +//_________________________________________________________________________________________________ +Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x) const +{ +// +// Implementation of the Background part of the Likelyhood function +// +// + Double_t retvalue = 0.; + retvalue = fParameters[7]*FunB(x) + (1. - fParameters[7])*FunP(x); + return retvalue; +} +//_________________________________________________________________________________________________ +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); + } + fitval = normFactorCB*fitvalCB; + return fitval; + }else{ + Double_t t=-1*m; + Double_t tmpv=-1*fParameters[9]; + fitval=TMath::Sqrt(TMath::Landau(t,tmpv,fParameters[11])); + fitval += fParameters[10]*(TMath::Landau(m,fParameters[9],fParameters[12])); + return fitval; + } +} +//_________________________________________________________________________________________________ +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.; + Double_t sigma3 = fResolutionConstants[2]; + Double_t xx; + Double_t sum = 0.0; + Double_t xlow,xupp; + Double_t step; + Double_t i; + 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); + + xx = xupp - (i-.5) * step; + sum += CsiMC(xx) * ResolutionFunc(xx-x); + } + return (step * sum) ; +} +//_________________________________________________________________________________________________ +Double_t AliBtoJPSItoEleCDFfitFCN::FunP(Double_t x) const +// +// Parameterisation of the Prompt part for the x distribution +// +{ + return ResolutionFunc(x); +} +//_________________________________________________________________________________________________ +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 +// + 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); +} +//_________________________________________________________________________________________________ +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); + 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); +} +//_________________________________________________________________________________________________ +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.; + Double_t sigma3 = fResolutionConstants[2]; + Double_t xx; + Double_t sum = 0.0; + Double_t xlow,xupp; + Double_t step; + Double_t i; + 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); + + 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); + } + 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.; + Double_t sigma3 = fResolutionConstants[2]; + Double_t xx; + Double_t sum = 0.0; + Double_t xlow,xupp; + Double_t step; + Double_t i; + 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); + } + 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.; + Double_t sigma3 = fResolutionConstants[2]; + Double_t xx; + Double_t sum1 = 0.0; + Double_t sum2 = 0.0; + Double_t xlow,xupp; + Double_t step; + Double_t i; + 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); + + 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); + } + 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); + } + 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]); +} diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h new file mode 100644 index 00000000000..297b61d4521 --- /dev/null +++ b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h @@ -0,0 +1,163 @@ +#ifndef ALIBTOJPSITOELECDFFITFCN_H +#define ALIBTOJPSITOELECDFFITFCN_H +/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//_________________________________________________________________________ +// Class AliBtoJPSItoEleCDFfitFCN +// Definition of main function used in +// unbinned log-likelihood fit for +// the channel B -> JPsi + X -> e+e- + X +// +// Origin: C.Di Giglio +// Contact: Carmelo.Digiglio@ba.infn.it , Giuseppe.Bruno@ba.infn.it +//_________________________________________________________________________ + +#include +#include +#include "TH1F.h" +#include "TMath.h" + +enum IntegralType {kBkg, + kBkgNorm, + kSig, + kSigNorm}; + +enum PtBins {kallpt, + kptbin1,kptbin2,kptbin3, + kptbin4,kptbin5,kptbin6, + kptbin7,kptbin8,kptbin9}; +//_________________________________________________________________________________________________ +class AliBtoJPSItoEleCDFfitFCN : public TNamed { + public: + // + AliBtoJPSItoEleCDFfitFCN(); + AliBtoJPSItoEleCDFfitFCN(const AliBtoJPSItoEleCDFfitFCN& source); + AliBtoJPSItoEleCDFfitFCN& operator=(const AliBtoJPSItoEleCDFfitFCN& source); + virtual ~AliBtoJPSItoEleCDFfitFCN(); + + 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;} + + void SetFPlus(Double_t plus) {fFPlus = plus;} + void SetFMinus(Double_t minus) {fFMinus = minus;} + void SetFSym(Double_t sym) {fFSym = sym;} + void SetRadius(Double_t radius) {fParameters[0] = radius;} + void SetTheta(Double_t theta) {fParameters[1] = theta;} + void SetPhi(Double_t phi) {fParameters[2] = phi;} + void SetLamPlus(Double_t lamplus) {fParameters[3] = lamplus;} + void SetLamMinus(Double_t lamminus) {fParameters[4] = lamminus;} + void SetLamSym(Double_t lamsym) {fParameters[5] = lamsym;} + void SetMassSlope(Double_t slope) {fParameters[6] = slope;} + void SetFractionJpsiFromBeauty(Double_t B) {fParameters[7] = B;} + void SetFsig(Double_t Fsig) {fParameters[8] = Fsig;} + void SetCrystalBallMmean(Double_t CrystalBallMmean) {fParameters[9] = CrystalBallMmean;} + 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 SetAllParameters(const Double_t* parameters); + void SetIntegral(Double_t integral) {fIntegral = 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 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.); + fFSym = TMath::Power((fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2])),2.);} + + void ComputeIntegral(); + + void ReadMCtemplates(Int_t BinNum); + + void PrintStatus(); + + private: + // + Double_t fParameters[13]; /* par[0] = fRadius; + par[1] = fTheta; + par[2] = fPhi; + par[3] = fOneOvLamPlus; + par[4] = fOneOvLamMinus; + par[5] = fOneOvLamSym; + par[6] = fMassBkgSlope; + par[7] = fFractionJpsiFromBeauty; + par[8] = fFsig; + 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 + + //// + + Double_t EvaluateCDFfunc(Double_t x, Double_t m) const ; + + Double_t EvaluateCDFfuncNorm(Double_t x, Double_t m) const ; + + //// + + Double_t EvaluateCDFfuncSignalPart(Double_t x, Double_t m) const ; // Signal part + + Double_t EvaluateCDFDecayTimeSigDistr(Double_t x) const ; + + Double_t EvaluateCDFInvMassSigDistr(Double_t m) const ; + + Double_t EvaluateCDFfuncBkgPart(Double_t x,Double_t m) const ; // Background part + + Double_t EvaluateCDFDecayTimeBkgDistr(Double_t x) const ; + + Double_t EvaluateCDFInvMassBkgDistr(Double_t m) const ; + + //// + + Double_t FunB(Double_t x) const ; + + Double_t FunP(Double_t x) const ; + + Double_t CsiMC(Double_t x) const ; + + Double_t FunBkgPos(Double_t x) const ; + + Double_t FunBkgNeg(Double_t x) const ; + + Double_t FunBkgSym(Double_t x) const ; + + Double_t ResolutionFunc(Double_t x) const ; + // + ClassDef (AliBtoJPSItoEleCDFfitFCN,0); // Unbinned log-likelihood fit + +}; + +#endif diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx new file mode 100644 index 00000000000..431e1509152 --- /dev/null +++ b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx @@ -0,0 +1,195 @@ +/************************************************************************** + * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ +#include +#include +#include "AliLog.h" +#include "AliBtoJPSItoEleCDFfitHandler.h" +#include "AliBtoJPSItoEleCDFfitFCN.h" + +//------------------------------------------------------------------------- +// Class AliBtoJPSItoEleCDFfitHandler +// Class to perform unbinned log-likelihood fit +// +// Origin: C. Di Giglio +// Contact: carmelo.digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it +//------------------------------------------------------------------------- + +void CDFFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag); + +ClassImp(AliBtoJPSItoEleCDFfitHandler) + +//______________________________________________________________________________ +void CDFFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) +{ + // This function is called by minuit + // The corresponding member method is called + // using SetObjectFit/GetObjectFit methods of TMinuit + AliBtoJPSItoEleCDFfitHandler* dummy = (AliBtoJPSItoEleCDFfitHandler *)TVirtualFitter::GetFitter()->GetObjectFit(); + dummy->CdfFCN(npar, gin, f, par, iflag); +} + + +//_________________________________________________________________________________________________ +AliBtoJPSItoEleCDFfitHandler::AliBtoJPSItoEleCDFfitHandler(): +fUp(0), +fX(0x0), +fM(0x0), +fLikely(0x0), +fNcand(0) +{ + // + // default constructor + // +} +//_________________________________________________________________________________________________ +AliBtoJPSItoEleCDFfitHandler::AliBtoJPSItoEleCDFfitHandler(Double_t* decaytime, + Double_t* invariantmass, Int_t ncand) : +fUp(0), +fX(decaytime), +fM(invariantmass), +fLikely(0x0), +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(Form("\n+++\n+++ Number of candidates ---> %d\n+++\n ", ncand)); +} +//___________________________________________________________________________ +AliBtoJPSItoEleCDFfitHandler& AliBtoJPSItoEleCDFfitHandler::operator=(const AliBtoJPSItoEleCDFfitHandler& c) +{ + // + // Assignment operator + // + if (this!=&c) { + fUp = c.fUp; + fX = c.fX; + fM = c.fM; + fLikely = c.fLikely; + fNcand = c.fNcand; + } + return *this; +} + +//___________________________________________________________________________ +AliBtoJPSItoEleCDFfitHandler::AliBtoJPSItoEleCDFfitHandler(const AliBtoJPSItoEleCDFfitHandler& c) : +TNamed(c), +fUp(c.fUp), +fX(c.fX), +fM(c.fM), +fLikely(c.fLikely), +fNcand(c.fNcand) +{ + // + // Copy Constructor + // +} +//________________________________________________________________________ +AliBtoJPSItoEleCDFfitHandler::~AliBtoJPSItoEleCDFfitHandler() +{ + // + //destructor + // + delete fLikely; +} +//_________________________________________________________________________________________________ +Int_t AliBtoJPSItoEleCDFfitHandler::DoMinimization() +{ + // + // performs the minimization + // + static TVirtualFitter *fitter = TVirtualFitter::Fitter(this,13); + 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}; + Int_t iret=fitter->ExecuteCommand("MIGRAD", arglist ,2); + fitter->PrintResults(4,0); + + AliInfo("Minimization procedure finished\n"); + + return iret; +} +//_______________________________________________________________________________________ +void AliBtoJPSItoEleCDFfitHandler::CdfFCN(Int_t & /* npar */, + Double_t * /* gin */,Double_t &f,Double_t *par,Int_t /* iflag */) +{ +// +// Definition of the FCN to be used by minuit +// + fLikely->SetAllParameters(par); + fLikely->ConvertFromSpherical(); + fLikely->ComputeIntegral(); + //printf("\n+++\n+++\n+++\n"); + //fLikely->PrintStatus(); + + TStopwatch t; + t.Start(); + + f = fLikely->EvaluateLikelihood(fX,fM,fNcand); + + t.Stop(); + AliDebug(2,Form("Real time spent to calculate function == %f \n", t.RealTime())); + AliDebug(2,Form("CPU time spent to calculate function == %f \n", t.CpuTime())); + AliDebug(2,Form("Actual value of the AliBtoJPSItoEleCDFfitFCN function == %f \n", f)); + + return; +} diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h new file mode 100644 index 00000000000..534934be30f --- /dev/null +++ b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h @@ -0,0 +1,49 @@ +#ifndef ALIBTOJPSITOELECDFFITHANDLER_H +#define ALIBTOJPSITOELECDFFITHANDLER_H +/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//___________________________________________________________________________ +// Class AliBtoJPSItoEleCDFfitHandler +// Class to perform unbinned log-likelihood fit +// +// Origin: C. Di Giglio +// Contact: Carmelo.Digiglio@ba.infn.it; Giuseppe.Bruno@ba.infn.it +//____________________________________________________________________________ + +#include +#include "TMath.h" +#include "AliBtoJPSItoEleCDFfitFCN.h" +//#include "AliLog.h" + +class AliBtoJPSItoEleCDFfitHandler : public TNamed { + public: + // + AliBtoJPSItoEleCDFfitHandler(); + AliBtoJPSItoEleCDFfitHandler& operator= (const AliBtoJPSItoEleCDFfitHandler& c); + 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()); } + void SetErrorDef(Double_t up) {fUp = up;} + + 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 */); + + Double_t* Decaytime() const { return fX; } + Double_t* InvariantMass() const { return fM; } + AliBtoJPSItoEleCDFfitFCN* LikelihoodPointer() const { return fLikely; } + Int_t DoMinimization(); + + private: + // + 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); + +}; +#endif -- 2.43.0