--- /dev/null
+/* 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;
+}
+//_________________________________________________________________________________________________
--- /dev/null
+#ifndef ALIANALYSISBTOJPSITOELE_H\r
+#define ALIANALYSISBTOJPSITOELE_H\r
+/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice */\r
+\r
+//-------------------------------------------------------------------------\r
+// Class AliAnalysisBtoJPSItoEle\r
+// Unbinned log-likelihood fit analysis class\r
+//\r
+// Origin: C.Di Giglio\r
+// Contact: Carmelo.Digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it\r
+//-------------------------------------------------------------------------\r
+class TNtuple ;\r
+class AliBtoJPSItoEleCDFfitFCN ;\r
+class AliBtoJPSItoEleCDFfitHandler ; \r
+#include "TH1F.h"\r
+\r
+//_________________________________________________________________________ \r
+class AliAnalysisBtoJPSItoEle : public TNamed {\r
+ public:\r
+ //\r
+ AliAnalysisBtoJPSItoEle();\r
+ AliAnalysisBtoJPSItoEle(const AliAnalysisBtoJPSItoEle& source);\r
+ AliAnalysisBtoJPSItoEle& operator=(const AliAnalysisBtoJPSItoEle& source);\r
+ virtual ~AliAnalysisBtoJPSItoEle();\r
+ \r
+ Int_t DoMinimization(Double_t* x,Double_t* m, Int_t n);\r
+ void ReadCandidates(TNtuple* nt, Double_t* &x, Double_t* &m, Int_t& n); // primary JPSI + secondary JPSI + bkg couples\r
+\r
+ void SetPtBin(Int_t BinNum) { fPtBin = BinNum ; }\r
+ void SetCsiMC(TH1F* MCtemplate);\r
+ void SetResolutionConstants(Int_t BinNum);\r
+ void CloneMCtemplate(const TH1F* MCtemplate) {fMCtemplate = (TH1F*)MCtemplate->Clone("fMCtemplate");}\r
+ Double_t* GetResolutionConstants(Double_t* resolutionConst);\r
+ Int_t GetPtBin() const { return fPtBin ; }\r
+\r
+ private:\r
+ //\r
+ AliBtoJPSItoEleCDFfitHandler* fFCNfunction; //! pointer to the interface class\r
+ Int_t fPtBin; // number of pt bin in which the analysis is performes\r
+ TH1F* fMCtemplate; //! template of the MC distribution for the x distribution of the secondary J/psi\r
+\r
+ ClassDef(AliAnalysisBtoJPSItoEle,0); // AliAnalysisBtoJPSItoEle class\r
+};\r
+\r
+#endif\r
--- /dev/null
+/**************************************************************************
+ * 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 <TFile.h>
+#include <TClonesArray.h>
+#include <TDatabasePDG.h>
+
+#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<AliAODEvent*> (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<TList*> (GetOutputData(1));
+ if (!fOutput) {
+ printf("ERROR: fOutput not available\n");
+ return;
+ }
+
+ if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
+ fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
+ fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
+ fhD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0"));
+ fhD0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0D0"));
+ fhCosThetaStar = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaStar"));
+ fhCosThetaPointing = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaPointing"));
+ fhCtsVsD0D0 = dynamic_cast<TH2F*>(fOutput->FindObject("fhCtsVsD0D0"));
+ fNtupleJPSI = dynamic_cast<TNtuple*>(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;
+
+}
+
--- /dev/null
+#ifndef ALIANALYSISTASKSEBTOJPSITOELE_H\r
+#define ALIANALYSISTASKSEBTOJPSITOELE_H\r
+/* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice */\r
+\r
+//*************************************************************************\r
+// Class AliAnalysisTaskSEBtoJPSItoEle\r
+// AliAnalysisTaskSE for the reconstruction \r
+// of heavy-flavour decay candidates\r
+// Author: C.Di Giglio, carmelo.digiglio@ba.infn.it\r
+//*************************************************************************\r
+\r
+class TH1F;\r
+class TList;\r
+class AliAnalysisBtoJPSItoEle;\r
+#include <TClonesArray.h>\r
+#include <TNtuple.h>\r
+#include <TH2F.h> \r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class AliAnalysisTaskSEBtoJPSItoEle : public AliAnalysisTaskSE\r
+{\r
+ public:\r
+\r
+ AliAnalysisTaskSEBtoJPSItoEle();\r
+ AliAnalysisTaskSEBtoJPSItoEle(const char *name);\r
+ virtual ~AliAnalysisTaskSEBtoJPSItoEle();\r
+\r
+ // Implementation of interface methods\r
+ virtual void UserCreateOutputObjects();\r
+ virtual void Init();\r
+ virtual void LocalInit() {Init();}\r
+ virtual void UserExec(Option_t *option);\r
+ virtual void Terminate(Option_t *option);\r
+\r
+ void SetCutsD0(const Double_t cutsD0[9]);\r
+ void SetPtCuts(const Double_t ptcuts[2]);\r
+ void SetAODMCInfo(Bool_t OkMCInfo) { fOkAODMC = OkMCInfo;}\r
+ void SetMinimize(Bool_t OkMinimize) { fOkMinimize = OkMinimize;}\r
+ void ReadAODMCInfo(const AliAODEvent* aodEv, const TClonesArray* inArray);\r
+ void SetPathMCfile (TString mcfilename="CsiMCfromKine_10PtBins.root") {fNameMCfile = mcfilename;}\r
+ TH1F* OpenLocalMCFile(TString dataDir, Int_t nbin); \r
+\r
+ private:\r
+\r
+ AliAnalysisTaskSEBtoJPSItoEle(const AliAnalysisTaskSEBtoJPSItoEle &source);\r
+ AliAnalysisTaskSEBtoJPSItoEle& operator=(const AliAnalysisTaskSEBtoJPSItoEle& source); \r
+ //\r
+ TList *fOutput; //! list send on output slot 0\r
+ AliAnalysisBtoJPSItoEle *fCdfFit; // Unbinned log-likelihood minimizer \r
+ TNtuple *fNtupleJPSI; // Ntuple of pseudo-proper decay time & invariant mass values\r
+ TH1F *fhDecayTimeMCjpsifromB; // Pseudo-proper decay time distribution used as template for JPSIs from B\r
+ TH1F *fhDecayTime; // Pseudo-proper decay time distribution\r
+ TH1F *fhInvMass; // Invariant mass distribution\r
+ TH1F *fhD0; // Impact parameter distribution\r
+ TH1F *fhD0D0; // Product of impact parameters distributions\r
+ TH1F *fhCosThetaStar; // Cosine of decay angle distribution\r
+ TH1F *fhCosThetaPointing; // Cosine of pointing angle distribution\r
+ TH2F *fhCtsVsD0D0; // Cos theta star Vs. D0D0 distribution\r
+ Bool_t fOkAODMC; // Flag to read AOD monte carlo information\r
+ TString fNameMCfile; // Name of the MC file with X template\r
+ Bool_t fOkMinimize; // Flag to enable unbinned log-likelihood minimization\r
+ \r
+ Double_t fCuts[9]; // cuts for N-tuple values selection\r
+ Double_t fPtCuts[2]; // Max and min pt of the candidates\r
+\r
+ ClassDef(AliAnalysisTaskSEBtoJPSItoEle,0); // AliAnalysisTaskSE for the reconstruction of heavy-flavour decay candidates\r
+};\r
+#endif\r
--- /dev/null
+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;
+
+}
--- /dev/null
+/**************************************************************************
+ * 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]);
+}
--- /dev/null
+#ifndef ALIBTOJPSITOELECDFFITFCN_H\r
+#define ALIBTOJPSITOELECDFFITFCN_H\r
+/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice */\r
+\r
+//_________________________________________________________________________\r
+// Class AliBtoJPSItoEleCDFfitFCN\r
+// Definition of main function used in \r
+// unbinned log-likelihood fit for\r
+// the channel B -> JPsi + X -> e+e- + X\r
+// \r
+// Origin: C.Di Giglio\r
+// Contact: Carmelo.Digiglio@ba.infn.it , Giuseppe.Bruno@ba.infn.it\r
+//_________________________________________________________________________\r
+\r
+#include <TNamed.h>\r
+#include <TDatabasePDG.h>\r
+#include "TH1F.h"\r
+#include "TMath.h"\r
+\r
+enum IntegralType {kBkg, \r
+ kBkgNorm, \r
+ kSig, \r
+ kSigNorm};\r
+\r
+enum PtBins {kallpt, \r
+ kptbin1,kptbin2,kptbin3,\r
+ kptbin4,kptbin5,kptbin6,\r
+ kptbin7,kptbin8,kptbin9};\r
+//_________________________________________________________________________________________________\r
+class AliBtoJPSItoEleCDFfitFCN : public TNamed {\r
+ public:\r
+ //\r
+ AliBtoJPSItoEleCDFfitFCN();\r
+ AliBtoJPSItoEleCDFfitFCN(const AliBtoJPSItoEleCDFfitFCN& source); \r
+ AliBtoJPSItoEleCDFfitFCN& operator=(const AliBtoJPSItoEleCDFfitFCN& source); \r
+ virtual ~AliBtoJPSItoEleCDFfitFCN();\r
+\r
+ Double_t EvaluateLikelihood(const Double_t* pseudoproperdecaytime,\r
+ const Double_t* invariantmass, const Int_t ncand);\r
+ \r
+ Double_t GetFPlus() const {return fFPlus;}\r
+ Double_t GetFMinus() const {return fFMinus;}\r
+ Double_t GetFSym() const {return fFSym;}\r
+ Double_t GetRadius() const {return fParameters[0];}\r
+ Double_t GetTheta() const {return fParameters[1];}\r
+ Double_t GetPhi() const {return fParameters[2];}\r
+ Double_t GetLamPlus() const {return fParameters[3];}\r
+ Double_t GetLamMinus() const {return fParameters[4];}\r
+ Double_t GetLamSym() const {return fParameters[5];}\r
+ Double_t GetMassSlope() const {return fParameters[6];}\r
+ Double_t GetFractionJpsiFromBeauty() const {return fParameters[7];}\r
+ Double_t GetFsig() const {return fParameters[8];}\r
+ Double_t GetCrystalBallMmean() const {return fParameters[9];}\r
+ Double_t GetCrystalBallNexp() const {return fParameters[10];}\r
+ Double_t GetCrystalBallSigma() const {return fParameters[11];}\r
+ Double_t GetCrystalBallAlpha() const {return fParameters[12];}\r
+ Double_t GetIntegral() const {return fIntegral;}\r
+ Bool_t GetCrystalBallParam() const {return fCrystalBallParam;}\r
+\r
+ void SetFPlus(Double_t plus) {fFPlus = plus;}\r
+ void SetFMinus(Double_t minus) {fFMinus = minus;}\r
+ void SetFSym(Double_t sym) {fFSym = sym;}\r
+ void SetRadius(Double_t radius) {fParameters[0] = radius;}\r
+ void SetTheta(Double_t theta) {fParameters[1] = theta;}\r
+ void SetPhi(Double_t phi) {fParameters[2] = phi;}\r
+ void SetLamPlus(Double_t lamplus) {fParameters[3] = lamplus;}\r
+ void SetLamMinus(Double_t lamminus) {fParameters[4] = lamminus;}\r
+ void SetLamSym(Double_t lamsym) {fParameters[5] = lamsym;}\r
+ void SetMassSlope(Double_t slope) {fParameters[6] = slope;}\r
+ void SetFractionJpsiFromBeauty(Double_t B) {fParameters[7] = B;}\r
+ void SetFsig(Double_t Fsig) {fParameters[8] = Fsig;}\r
+ void SetCrystalBallMmean(Double_t CrystalBallMmean) {fParameters[9] = CrystalBallMmean;}\r
+ void SetCrystalBallNexp(Double_t CrystalBallNexp) {fParameters[10] = CrystalBallNexp;}\r
+ void SetCrystalBallSigma(Double_t CrystalBallSigma) {fParameters[11] = CrystalBallSigma;}\r
+ void SetCrystalBallAlpha(Double_t CrystalBallAlpha) {fParameters[12] = CrystalBallAlpha;}\r
+\r
+ void SetAllParameters(const Double_t* parameters);\r
+ void SetIntegral(Double_t integral) {fIntegral = integral;}\r
+ void SetCsiMC(const TH1F* MCtemplate) {fhCsiMC = (TH1F*)MCtemplate->Clone("fhCsiMC");}\r
+ void SetResolutionConstants(Int_t BinNum);\r
+ void SetMassWndHigh(Double_t limit) { fMassWndHigh = TDatabasePDG::Instance()->GetParticle(443)->Mass() + limit ;}//here use pdg code instead\r
+ void SetMassWndLow(Double_t limit) { fMassWndLow = TDatabasePDG::Instance()->GetParticle(443)->Mass() - limit ;}//here use pdg code instead\r
+ void SetCrystalBallParam(Bool_t okCB) {fCrystalBallParam = okCB;}\r
+\r
+ void ConvertFromSpherical() { fFPlus = TMath::Power((fParameters[0]*TMath::Cos(fParameters[1])),2.);\r
+ fFMinus = TMath::Power((fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2])),2.);\r
+ fFSym = TMath::Power((fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2])),2.);} \r
+\r
+ void ComputeIntegral(); \r
+\r
+ void ReadMCtemplates(Int_t BinNum);\r
+\r
+ void PrintStatus();\r
+\r
+ private: \r
+ //\r
+ Double_t fParameters[13]; /* par[0] = fRadius; \r
+ par[1] = fTheta;\r
+ par[2] = fPhi;\r
+ par[3] = fOneOvLamPlus;\r
+ par[4] = fOneOvLamMinus;\r
+ par[5] = fOneOvLamSym;\r
+ par[6] = fMassBkgSlope;\r
+ par[7] = fFractionJpsiFromBeauty;\r
+ par[8] = fFsig;\r
+ par[9] = fCrystalBallMmean;\r
+ par[10] = fCrystalBallNexp;\r
+ par[11] = fCrystalBallSigma;\r
+ par[12] = fCrystalBallAlpha;*/\r
+\r
+ Double_t fFPlus; // parameters of the log-likelihood function\r
+ Double_t fFMinus; // Slopes of the x distributions of the background \r
+ Double_t fFSym; // functions \r
+\r
+ Double_t fIntegral; // integral values of log-likelihood function terms\r
+ TH1F *fhCsiMC; // X distribution used as MC template for JPSI from B\r
+ Double_t fResolutionConstants[7]; // constants for the parametrized resolution function R(X)\r
+ Double_t fMassWndHigh; // JPSI Mass window higher limit\r
+ Double_t fMassWndLow; // JPSI Mass window lower limit\r
+ Bool_t fCrystalBallParam; // Boolean to switch to Crystall Ball parameterisation\r
+\r
+ ////\r
+\r
+ Double_t EvaluateCDFfunc(Double_t x, Double_t m) const ;\r
+\r
+ Double_t EvaluateCDFfuncNorm(Double_t x, Double_t m) const ;\r
+\r
+ ////\r
+\r
+ Double_t EvaluateCDFfuncSignalPart(Double_t x, Double_t m) const ; // Signal part \r
+\r
+ Double_t EvaluateCDFDecayTimeSigDistr(Double_t x) const ;\r
+\r
+ Double_t EvaluateCDFInvMassSigDistr(Double_t m) const ;\r
+\r
+ Double_t EvaluateCDFfuncBkgPart(Double_t x,Double_t m) const ; // Background part\r
+\r
+ Double_t EvaluateCDFDecayTimeBkgDistr(Double_t x) const ;\r
+ \r
+ Double_t EvaluateCDFInvMassBkgDistr(Double_t m) const ;\r
+\r
+ ////\r
+\r
+ Double_t FunB(Double_t x) const ;\r
+\r
+ Double_t FunP(Double_t x) const ;\r
+\r
+ Double_t CsiMC(Double_t x) const ;\r
+\r
+ Double_t FunBkgPos(Double_t x) const ;\r
+\r
+ Double_t FunBkgNeg(Double_t x) const ;\r
+\r
+ Double_t FunBkgSym(Double_t x) const ;\r
+\r
+ Double_t ResolutionFunc(Double_t x) const ;\r
+ //\r
+ ClassDef (AliBtoJPSItoEleCDFfitFCN,0); // Unbinned log-likelihood fit \r
+\r
+};\r
+\r
+#endif\r
--- /dev/null
+/**************************************************************************
+ * 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 <TVirtualFitter.h>
+#include <TStopwatch.h>
+#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;
+}
--- /dev/null
+#ifndef ALIBTOJPSITOELECDFFITHANDLER_H\r
+#define ALIBTOJPSITOELECDFFITHANDLER_H\r
+/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice */\r
+\r
+//___________________________________________________________________________\r
+// Class AliBtoJPSItoEleCDFfitHandler\r
+// Class to perform unbinned log-likelihood fit\r
+// \r
+// Origin: C. Di Giglio\r
+// Contact: Carmelo.Digiglio@ba.infn.it; Giuseppe.Bruno@ba.infn.it\r
+//____________________________________________________________________________\r
+\r
+#include <TNamed.h>\r
+#include "TMath.h"\r
+#include "AliBtoJPSItoEleCDFfitFCN.h"\r
+//#include "AliLog.h"\r
+\r
+class AliBtoJPSItoEleCDFfitHandler : public TNamed {\r
+ public:\r
+ //\r
+ AliBtoJPSItoEleCDFfitHandler();\r
+ AliBtoJPSItoEleCDFfitHandler& operator= (const AliBtoJPSItoEleCDFfitHandler& c);\r
+ AliBtoJPSItoEleCDFfitHandler(const AliBtoJPSItoEleCDFfitHandler& c);\r
+ AliBtoJPSItoEleCDFfitHandler(Double_t* decaytime, Double_t* invariantmass, Int_t ncand);\r
+ ~AliBtoJPSItoEleCDFfitHandler(); \r
+ Double_t Up() const { return fUp*TMath::Sqrt(fLikely->GetIntegral()); }\r
+ void SetErrorDef(Double_t up) {fUp = up;}\r
+\r
+ Double_t operator()(const Double_t* par) const ;\r
+ void CdfFCN(Int_t & /* npar */, Double_t * /* gin */, Double_t &f, Double_t *par, Int_t /* iflag */);\r
+\r
+ Double_t* Decaytime() const { return fX; }\r
+ Double_t* InvariantMass() const { return fM; }\r
+ AliBtoJPSItoEleCDFfitFCN* LikelihoodPointer() const { return fLikely; }\r
+ Int_t DoMinimization();\r
+\r
+ private:\r
+ //\r
+ Double_t fUp; //error definition \r
+ Double_t* fX; //pseudo-proper decay time X\r
+ Double_t* fM; //invariant mass M\r
+ AliBtoJPSItoEleCDFfitFCN* fLikely; //Log likelihood function\r
+ Int_t fNcand; //number of candidates\r
+ //\r
+ ClassDef(AliBtoJPSItoEleCDFfitHandler,0);\r
+\r
+}; \r
+#endif\r
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff