o updates for PbPb analysis of B->J/psi (Fiorella)
authorwiechula <wiechula@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 3 Nov 2013 22:53:46 +0000 (22:53 +0000)
committerwiechula <wiechula@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 3 Nov 2013 22:53:46 +0000 (22:53 +0000)
PWGDQ/dielectron/AliAnalysisTaskDielectronReadAODBranch.cxx
PWGDQ/dielectron/AliAnalysisTaskDielectronReadAODBranch.h
PWGDQ/dielectron/AliDielectronBtoJPSItoEle.cxx
PWGDQ/dielectron/AliDielectronBtoJPSItoEle.h
PWGDQ/dielectron/AliDielectronBtoJPSItoEleCDFfitFCN.cxx
PWGDQ/dielectron/AliDielectronBtoJPSItoEleCDFfitFCN.h
PWGDQ/dielectron/AliDielectronBtoJPSItoEleCDFfitFCNfitter.cxx
PWGDQ/dielectron/AliDielectronBtoJPSItoEleCDFfitHandler.cxx
PWGDQ/dielectron/AliDielectronBtoJPSItoEleCDFfitHandler.h

index b5dc260..7cf300a 100644 (file)
@@ -22,6 +22,7 @@
 
 #include "TChain.h"
 #include "TNtuple.h"
+#include "TList.h"
 #include "TH1F.h"
 #include "TH2F.h"
 #include "TDatabasePDG.h"
@@ -34,6 +35,7 @@
 #include "AliAODMCParticle.h"
 #include "AliAnalysisTaskDielectronReadAODBranch.h"
 #include "AliDielectronPair.h"
+#include "AliDielectronMC.h"
 
 ClassImp(AliAnalysisTaskDielectronReadAODBranch)
        //________________________________________________________________________
@@ -73,7 +75,8 @@ ClassImp(AliAnalysisTaskDielectronReadAODBranch)
                 fPairType(1),
                 fPtJpsi(1.3),
                 fInvMassSignalLimits(0),
-                fInvMassSideBandsLimits(0)
+                fInvMassSideBandsLimits(0),
+               fSecondary(0)
 {
        // Default constructor
 }
@@ -115,7 +118,8 @@ AliAnalysisTaskDielectronReadAODBranch::AliAnalysisTaskDielectronReadAODBranch(c
         fPairType(1),
         fPtJpsi(1.3),
         fInvMassSignalLimits(0),
-        fInvMassSideBandsLimits(0) 
+        fInvMassSideBandsLimits(0), 
+       fSecondary(0)
 {
        // Default constructor
        DefineInput(0,TChain::Class());
@@ -158,7 +162,7 @@ void AliAnalysisTaskDielectronReadAODBranch::UserCreateOutputObjects()
         fInvMass = new TH1F("fInvMass","fInvMass",300,2.0,2.0+300*.04); // step 40MeV
         fInvMassNoCuts = new TH1F("fInvMass_no_cuts","fInvMass_no_cuts",125,0,125*0.04); //step 40MeV
         fNentries = new TH1F("numberOfevent","numbOfEvent",1,-0.5,0.5);
-       
+
         //pseudoproper 
        fpsproperSignal = new TH1F("psproper_decay_length",Form("psproper_decay_length_distrib(AODcuts+#clustTPC>%d, pT(leg)>%fGeV, pT(J/#psi)>%fGeV/c, %f < M < %f);X [#mu m];Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]),150,-3000.,3000.);
        fpsproperSidebands = new TH1F("psproper_decay_length_sidebands",Form("psproper_decay_length_distrib_sidebands(AODcuts+#clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,M> %f && M < %f );  X [#mu m]; Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSideBandsLimits[1],fInvMassSideBandsLimits[0]),150,-3000.,3000.);
@@ -223,7 +227,7 @@ void AliAnalysisTaskDielectronReadAODBranch::UserExec(Option_t */*option*/)
        // Execute analysis for current event:
        AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
        if (!aod) return;
-       
+
        Double_t vtxPrim[3] = {0.,0.,0.};
 
        AliAODVertex* primvtx = aod->GetPrimaryVertex();
@@ -241,8 +245,6 @@ void AliAnalysisTaskDielectronReadAODBranch::UserExec(Option_t */*option*/)
 
        fNentries->Fill(0);
        aodTree->GetEvent(Entry());
-       Int_t dgLabels[2] = {0,0};
-       Int_t pdgDgJPSItoEE[2] = {11,11};
 
        // loop over candidates
        if(fobj) {
@@ -266,12 +268,10 @@ void AliAnalysisTaskDielectronReadAODBranch::UserExec(Option_t */*option*/)
                        AliAODTrack *trk = (AliAODTrack*)pairObj->GetFirstDaughter();
                        AliAODTrack *trk1 = (AliAODTrack*)pairObj->GetSecondDaughter();
                        if(!trk || !trk1) {printf("ERROR: daughter tracks not available\n"); continue;}
-                       dgLabels[0] = trk->GetLabel(); 
-                       dgLabels[1] = trk1->GetLabel(); 
 
                        //check in case of MC analysis if candidate is a true J/psi->ee
-                       if(fHasMC && !MatchToMC(443,fobjMC,dgLabels,2,2,pdgDgJPSItoEE)) continue;
-                           AliAODPid *pid  = trk->GetDetPid();
+                          if(fHasMC && ((AliDielectronMC::Instance()->IsJpsiPrimary(pairObj))!=fSecondary)) continue;
+                          AliAODPid *pid  = trk->GetDetPid();
                                 AliAODPid *pid1 = trk1->GetDetPid();
                                 if(!pid || !pid1){
                                         if(fDebug>1) printf("No AliAODPid found\n");
@@ -337,122 +337,3 @@ void AliAnalysisTaskDielectronReadAODBranch::Terminate(Option_t */*option*/)
        return;
 }
 
-//___________________________________________________________________
-Bool_t AliAnalysisTaskDielectronReadAODBranch::MatchToMC(Int_t pdgabs,const TObjArray *mcArray,
-               const Int_t *dgLabels,Int_t ndg,
-               Int_t ndgCk, const Int_t *pdgDg) const
-{
-       //
-       // jpsi candidates association to MC truth
-        // Check if this candidate is matched to a MC signal
-       // If no, return kFALSE
-       // If yes, return kTRUE
-       //
-       if (!mcArray) return kFALSE;
-       
-       Int_t labMom[2]={0,0};
-       Int_t i,j,lab,labMother,pdgMother,pdgPart,labJPSIMother,pdgJPSIMother;
-       AliAODMCParticle *part=0;
-       AliAODMCParticle *mother=0;
-       AliAODMCParticle *jPSImother=0;
-       Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.;
-       Bool_t pdgUsed[2]={kFALSE,kFALSE};
-
-       // loop on daughter labels
-       for(i=0; i<ndg; i++) {
-               labMom[i]=-1;
-               lab = dgLabels[i];
-               if(lab<0) {
-                       printf("daughter with negative label %d\n",lab);
-                       return kFALSE;
-               }
-               part = (AliAODMCParticle*)mcArray->At(lab);
-               if(!part) {
-                       printf("no MC particle\n");
-                       return kFALSE;
-               }
-
-               // check the PDG of the daughter, if requested
-               if(ndgCk>0) {
-                       pdgPart=TMath::Abs(part->GetPdgCode());
-                       printf("pdg code of daughter %d == %d\n",i,pdgPart);
-                       for(j=0; j<ndg; j++) {
-                               if(!pdgUsed[j] && pdgPart==pdgDg[j]) {
-                                       pdgUsed[j]=kTRUE;
-                                       break;
-                               }
-                       }
-               }
-
-               // for the J/psi, check that the daughters are electrons
-               if(pdgabs==443) {
-                       if(TMath::Abs(part->GetPdgCode())!=11) return kFALSE;
-               }
-
-               mother = part;
-               while(mother->GetMother()>=0) {
-                       labMother=mother->GetMother();
-                       mother = (AliAODMCParticle*)mcArray->At(labMother);//get particle mother
-                       if(!mother) {
-                               printf("no MC mother particle\n");
-                               break;
-                       }
-                       pdgMother = TMath::Abs(mother->GetPdgCode());//get particle mother pdg code
-                       printf("pdg code of mother track == %d\n",pdgMother);
-                       if(pdgMother==pdgabs) {
-                               labJPSIMother=mother->GetMother();
-                               jPSImother = (AliAODMCParticle*)mcArray->At(labJPSIMother);//get J/psi mother
-                               if(jPSImother) {
-                                       pdgJPSIMother = TMath::Abs(jPSImother->GetPdgCode()); //get J/psi mother pdg code
-                                       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) return kFALSE;} //check if J/psi comes from B hadron
-
-                                       labMom[i]=labMother;
-                                       // keep sum of daughters' momenta, to check for mom conservation
-                                       pxSumDgs += part->Px();
-                                       pySumDgs += part->Py();
-                                       pzSumDgs += part->Pz();
-                                       break;
-                       } else if(pdgMother>pdgabs || pdgMother<10) {
-                               break;
-                       }
-               }
-               if(labMom[i]==-1) {printf("mother PDG not ok for this daughter\n"); return kFALSE;} // mother PDG not ok for this daughter
-       } // end loop on daughters
-
-       // check if the candidate is signal
-       labMother=labMom[0];
-       // all labels have to be the same and !=-1
-       for(i=0; i<ndg; i++) {
-               if(labMom[i]==-1)        return kFALSE;
-               if(labMom[i]!=labMother) return kFALSE;
-       }
-
-       // check that all daughter PDGs are matched
-       if(ndgCk>0) {
-               for(i=0; i<ndg; i++) {
-                       if(pdgUsed[i]==kFALSE) return kFALSE;
-               }
-       }
-
-       mother = (AliAODMCParticle*)mcArray->At(labMother);
-       Double_t pxMother = mother->Px();
-       Double_t pyMother = mother->Py();
-       Double_t pzMother = mother->Pz();
-       // within 0.1%
-       if((TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) > 0.00001 &&
-                       (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.00001 &&
-                       (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.00001) return kFALSE;
-
-       return kTRUE;
-}
index 6e712c9..a350b4e 100644 (file)
@@ -20,6 +20,7 @@
 class TNtuple;
 class TH1F;
 class TH2F;
+class AliDielectronPair;
 
 class AliAnalysisTaskDielectronReadAODBranch : public AliAnalysisTaskSE
 {
@@ -34,7 +35,6 @@ class AliAnalysisTaskDielectronReadAODBranch : public AliAnalysisTaskSE
                virtual void LocalInit() {Init();}
                virtual void UserExec(Option_t *option);
                virtual void Terminate(Option_t *option);
-               Bool_t MatchToMC(Int_t pdgabs,const TObjArray *mcArray,const Int_t *dgLabels,Int_t ndg,Int_t ndgCk,const Int_t *pdgDg) const;
               
                 //setters
                 void SetPtLeg(Double_t cutPt){ fPtCut = cutPt;} 
@@ -99,6 +99,7 @@ class AliAnalysisTaskDielectronReadAODBranch : public AliAnalysisTaskSE
                 Double_t  *fInvMassSignalLimits;     // invariant mass signal region to extract psproper distribution
                 Double_t  *fInvMassSideBandsLimits;  // invariant mass sideband region to extract psproper distribution 
 
+                Int_t     fSecondary;                // 1(0) to select secondary (prompt) jpsi
                ClassDef(AliAnalysisTaskDielectronReadAODBranch,2); // AliAnalysisTaskSE for the MC association of heavy-flavour decay candidates
 };
 
index 8550237..a02bdf7 100644 (file)
@@ -22,7 +22,6 @@ ClassImp(AliDielectronBtoJPSItoEle)
 //_______________________________________________________________________________ 
 AliDielectronBtoJPSItoEle::AliDielectronBtoJPSItoEle() :
 fFCNfunction(0),
-fPtBin(0),
 fMCtemplate(0),
 fResType("FF")
 {
@@ -34,7 +33,6 @@ fResType("FF")
 AliDielectronBtoJPSItoEle::AliDielectronBtoJPSItoEle(const AliDielectronBtoJPSItoEle& source) :
 TNamed(source),
 fFCNfunction(source.fFCNfunction),
-fPtBin(source.fPtBin),
 fMCtemplate(source.fMCtemplate),
 fResType(source.fResType)
 {
@@ -51,7 +49,6 @@ AliDielectronBtoJPSItoEle &AliDielectronBtoJPSItoEle::operator=(const AliDielect
   //
   if(&source == this) return *this;
   fFCNfunction = source.fFCNfunction;
-  fPtBin = source.fPtBin;
   fMCtemplate = source.fMCtemplate;
 
   return *this;
@@ -76,17 +73,18 @@ Int_t AliDielectronBtoJPSItoEle::DoMinimization(Int_t step)
   return iret;
 }
 //_________________________________________________________________________________________________
-void AliDielectronBtoJPSItoEle::ReadCandidates(TNtuple* nt, Double_t* &pseudoproper, Double_t* &invmass, Int_t * &typeCand, Int_t& ncand)
+void AliDielectronBtoJPSItoEle::ReadCandidates(TNtuple* nt, Double_t* &pseudoproper, Double_t* &invmass, Double_t* &pt, Int_t * &typeCand, Int_t& ncand, Double_t massLow, Double_t massUp, Double_t ptLow, Double_t ptUp)
 {
   //
   // Read N-tuple with X and M values
   //
-  Float_t mJPSI = 0; Float_t x = 0; Float_t type = 0;
+  Float_t mJPSI = 0; Float_t x = 0; Float_t type = 0; Float_t transvMom = 0.;
   Int_t nentries = 0;
   ncand=0;
   TString arrType[] = {"SS","FS","FF"};
   nt->SetBranchAddress("Mass",&mJPSI);
   nt->SetBranchAddress("Xdecaytime",&x);
+  nt->SetBranchAddress("Pt",&transvMom);
   //
   if(!nt->GetListOfBranches()->At(2)) {AliInfo("ERROR: branch with candidate type doesn't exist! \n"); return;}
   nt->SetBranchAddress("Type",&type);
@@ -94,13 +92,17 @@ void AliDielectronBtoJPSItoEle::ReadCandidates(TNtuple* nt, Double_t* &pseudopro
   nentries = (Int_t)nt->GetEntries();
   pseudoproper = new Double_t[nentries];
   invmass      = new Double_t[nentries];
+  pt           = new Double_t[nentries];
   typeCand     = new Int_t[nentries];
 
   for(Int_t i = 0; i < nentries; i++) {
       nt->GetEntry(i);
       if(!fResType.Contains(arrType[(Int_t)type])) continue;
+      if(massUp > massLow && massLow > 0) { if(mJPSI < massLow || mJPSI >massUp) continue; } 
+      if(ptUp > ptLow && ptLow > 0) { if(transvMom < ptLow || transvMom > ptUp) continue; } 
       pseudoproper[ncand]=(Double_t)x;
       invmass[ncand]=(Double_t)mJPSI;
+      pt[ncand]=(Double_t)transvMom;
       typeCand[ncand] = (Int_t)type;
       ncand++;
       }
@@ -117,13 +119,13 @@ void AliDielectronBtoJPSItoEle::SetCsiMC()
   return;
 }
 //_________________________________________________________________________________________________
-void AliDielectronBtoJPSItoEle::SetFitHandler(Double_t* x /*pseudoproper*/, Double_t* m /*inv mass*/, Int_t* type /*type*/, Int_t ncand /*candidates*/) 
+void AliDielectronBtoJPSItoEle::SetFitHandler(Double_t* x /*pseudoproper*/, Double_t* m /*inv mass*/, Double_t* pt /*pt*/, Int_t* type /*type*/, Int_t ncand /*candidates*/) 
 {
   //
   // Create the fit handler object to play with different params of the fitting function
   //
 
-  fFCNfunction = new AliDielectronBtoJPSItoEleCDFfitHandler(x,m,type,ncand);
+  fFCNfunction = new AliDielectronBtoJPSItoEleCDFfitHandler(x,m,pt,type,ncand);
   if(!fFCNfunction) {
 
      AliInfo("fFCNfunction not istanziated  ---> nothing done");
index e2ef570..016c9c2 100644 (file)
@@ -26,21 +26,18 @@ class AliDielectronBtoJPSItoEle : public TNamed {
                virtual ~AliDielectronBtoJPSItoEle();\r
 \r
                Int_t DoMinimization(Int_t step = 0);\r
-               void ReadCandidates(TNtuple* nt, Double_t* &x, Double_t* &m, Int_t * &typeCand, Int_t& n); // primary JPSI + secondary JPSI + bkg couples\r
+               void ReadCandidates(TNtuple* nt, Double_t* &x, Double_t* &m, Double_t* &pt, Int_t * &typeCand, Int_t& n,Double_t massLow = -1., Double_t massUp = -1., Double_t ptLow = -1., Double_t ptUp = -1.); // primary JPSI + secondary JPSI + bkg couples\r
 \r
-               void SetPtBin(Int_t BinNum) { fPtBin = BinNum ; }\r
                void SetCsiMC();\r
-               void SetFitHandler(Double_t* x /*pseudoproper*/, Double_t* m /*inv mass*/, Int_t *type /*type*/, Int_t ncand /*candidates*/); \r
+               void SetFitHandler(Double_t* x /*pseudoproper*/, Double_t* m /*inv mass*/, Double_t *pt /*transverse momentum */, Int_t *type /*type*/, Int_t ncand /*candidates*/); \r
                void CloneMCtemplate(const TH1F* MCtemplate) {fMCtemplate = (TH1F*)MCtemplate->Clone("fMCtemplate");}\r
                void SetResTypeAnalysis(TString resType){fResType = resType;}\r
                 Double_t* GetResolutionConstants(Double_t* resolutionConst);\r
                AliDielectronBtoJPSItoEleCDFfitHandler* GetCDFFitHandler() const { return fFCNfunction ; }\r
-               Int_t GetPtBin() const { return fPtBin ; }\r
 \r
        private:\r
                //\r
                AliDielectronBtoJPSItoEleCDFfitHandler* 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
                 TString fResType;                           // string with candidate's types considered\r
 \r
index 72ff24c..04aa8bc 100644 (file)
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-#include "AliLog.h"
+#include <AliLog.h>
+#include <TFormula.h>
+#include <TF1.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TFile.h>
+
 #include "AliDielectronBtoJPSItoEleCDFfitFCN.h"
-#include "TFormula.h"
-#include "TF1.h"
-#include "TCanvas.h"
-#include "TMath.h"
 
 //_________________________________________________________________________
 //                        Class AliDielectronBtoJPSItoEleCDFfitFCN
@@ -39,38 +41,67 @@ ClassImp(AliDielectronBtoJPSItoEleCDFfitFCN)
                fintmMassSig(1.),
                fintmMassBkg(1.),
                fhCsiMC(0x0),
+                fShiftTemplate(0.),
                fMassWndHigh(0.),
                fMassWndLow(0.),
-               fCrystalBallParam(kFALSE)
+               fCrystalBallParam(kFALSE),
+                fChangeResolution(1.),
+                fChangeMass(1.),
+               fWeights(0),
+                fLoadFunctions(kFALSE),
+                fMultivariate(kFALSE),
+                fFunBSaved(0x0),
+                fFunBkgSaved(0x0),
+                fResParams(0x0),
+                fBkgParams(0x0),
+               fMassWindows(0x0),
+               fPtWindows(0x0),
+               fExponentialParam(kTRUE),
+               fSignalBinForExtrapolation(0)
 {
        //
        // constructor
        //
-       SetCrystalBallFunction(kFALSE);
        SetMassWndHigh(0.2);
        SetMassWndLow(0.5);
         fWeightType[0] = 1.; fWeightType[1] = 1.; fWeightType[2] = 1.;
-        for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = 0.;
+        for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = 0.;
         fParameters[9] = 1.;fParameters[11] = 1.;fParameters[12] = 1.;
-       AliInfo("Instance of AliDielectronBtoJPSItoEleCDFfitFCN-class created");
+        
+        
+        AliInfo("Instance of AliDielectronBtoJPSItoEleCDFfitFCN-class created");
 }
 //_________________________________________________________________________________________________
 AliDielectronBtoJPSItoEleCDFfitFCN::AliDielectronBtoJPSItoEleCDFfitFCN(const AliDielectronBtoJPSItoEleCDFfitFCN& source) :
        TNamed(source),
        fFPlus(source.fFPlus),
-       fFMinus(source.fFMinus),
-       fFSym(source.fFSym),
-       fintmMassSig(source.fintmMassSig),
-       fintmMassBkg(source.fintmMassBkg),
-       fhCsiMC(source.fhCsiMC),
-       fMassWndHigh(source.fMassWndHigh),
-       fMassWndLow(source.fMassWndLow),
-       fCrystalBallParam(source.fCrystalBallParam) 
+        fFMinus(source.fFMinus),
+        fFSym(source.fFSym),
+        fintmMassSig(source.fintmMassSig),
+        fintmMassBkg(source.fintmMassBkg),
+        fhCsiMC(source.fhCsiMC),
+        fShiftTemplate(source.fShiftTemplate),
+        fMassWndHigh(source.fMassWndHigh),
+        fMassWndLow(source.fMassWndLow),
+        fCrystalBallParam(source.fCrystalBallParam),
+        fChangeResolution(source.fChangeResolution),
+        fChangeMass(source.fChangeMass),
+        fWeights(source.fWeights),
+        fLoadFunctions(source.fLoadFunctions),
+        fMultivariate(source.fMultivariate),
+        fFunBSaved(source.fFunBSaved),
+        fFunBkgSaved(source.fFunBkgSaved),
+        fResParams(source.fResParams),
+        fBkgParams(source.fBkgParams),
+        fMassWindows(source.fMassWindows),
+        fPtWindows(source.fPtWindows),
+        fExponentialParam(source.fExponentialParam),
+        fSignalBinForExtrapolation(source.fSignalBinForExtrapolation)
 {
        //
        // Copy constructor
        //
-       for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = source.fParameters[iPar];
+       for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = source.fParameters[iPar];
         for(Int_t iW=0; iW<2; iW++) fWeightType[iW] = source.fWeightType[iW]; 
 }
 //_________________________________________________________________________________________________
@@ -80,15 +111,25 @@ AliDielectronBtoJPSItoEleCDFfitFCN& AliDielectronBtoJPSItoEleCDFfitFCN::operator
        // Assignment operator
        //
        if(&source == this) return *this;
-       fFPlus = source.fFPlus;
+       fFPlus = source.fFPlus;
        fFMinus = source.fFMinus;
        fFSym = source.fFSym;
        fintmMassSig = source.fintmMassSig;
        fintmMassBkg = source.fintmMassBkg;
        fhCsiMC = source.fhCsiMC;
+        fLoadFunctions = source.fLoadFunctions;
+        fMultivariate = source.fMultivariate;
+        fFunBSaved = source.fFunBSaved;
+        fFunBkgSaved = source.fFunBkgSaved;
+        fResParams = source.fResParams;
+        fBkgParams = source.fBkgParams;
+        fMassWindows = source.fMassWindows;
+        fPtWindows = source.fPtWindows;
+        fShiftTemplate = source.fShiftTemplate;
        fCrystalBallParam = source.fCrystalBallParam;
+       fExponentialParam = source.fExponentialParam;
 
-       for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = source.fParameters[iPar];
+       for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = source.fParameters[iPar];
        return *this;
 }  
 //_________________________________________________________________________________________________
@@ -99,11 +140,13 @@ AliDielectronBtoJPSItoEleCDFfitFCN::~AliDielectronBtoJPSItoEleCDFfitFCN()
        //
 
        delete fhCsiMC;
-       for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = 0.;
+        delete fFunBSaved;
+        delete fFunBkgSaved;
+       for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = 0.;
 }
 //_________________________________________________________________________________________________
 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t* pseudoproperdecaytime,
-               const Double_t* invariantmass, const Int_t *type, const Int_t ncand) const
+               const Double_t* invariantmass, const Double_t *pt, const Int_t *type, const Int_t ncand) const
 {
        //
        // This function evaluates the Likelihood fnction
@@ -113,9 +156,9 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t*
        Double_t ret = 0.;
 
        for(Int_t i=0; i < ncand; i++) {
-               f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i],type[i]);
-               if(f <= 0.) continue;  
-                ret += -1*TMath::Log(f);  
+               f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i],pt[i],type[i]);
+               if(f <= 0.) continue;   
+                ret += -2.*TMath::Log(f);  
        }
         return ret;
 }
@@ -125,7 +168,7 @@ void AliDielectronBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parame
        //
        // Sets array of FCN parameters
        //
-       for(Int_t index = 0; index < 45; index++) fParameters[index] = parameters[index];
+       for(Int_t index = 0; index < 49; index++) fParameters[index] = parameters[index];
 }
 //_________________________________________________________________________________________________
 void AliDielectronBtoJPSItoEleCDFfitFCN::ComputeMassIntegral() 
@@ -177,9 +220,11 @@ void AliDielectronBtoJPSItoEleCDFfitFCN::PrintStatus()
        printf("actual value of fPos ------------------------------------------>> | %f \n", GetFPlus());
        printf("actual value of fNeg ------------------------------------------>> | %f \n", GetFMinus());
        printf("actual value of fSym ------------------------------------------>> | %f \n", GetFSym()); 
+       printf("actual value of fSym1 ----------------------------------------->> | %f \n", GetFSym1()); 
        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 fOneOvLamSym1 --------------------------------->> | %f \n", GetLamSym1());
        printf("actual value of fFractionJpsiFromBeauty ----------------------->> | %f \n", GetFractionJpsiFromBeauty());
        printf("actual value of fFsig ----------------------------------------->> | %f \n", GetFsig());
 
@@ -197,10 +242,20 @@ void AliDielectronBtoJPSItoEleCDFfitFCN::PrintStatus()
        }
 
        // back Mass func
-       printf("actual value of normBkg ----------------------------------------->> | %f \n", GetBkgInvMassNorm());
-       printf("actual value of meanBkg ----------------------------------------->> | %f \n", GetBkgInvMassMean());
-       printf("actual value of slopeBkg ---------------------------------------->> | %f \n", GetBkgInvMassSlope());
-       printf("actual value of constBkg ---------------------------------------->> | %f \n", GetBkgInvMassConst());
+        if(fExponentialParam){
+               printf("actual value of normBkg ----------------------------------------->> | %f \n", GetBkgInvMassNorm());
+               printf("actual value of meanBkg ----------------------------------------->> | %f \n", GetBkgInvMassMean());
+               printf("actual value of slopeBkg ---------------------------------------->> | %f \n", GetBkgInvMassSlope());
+               printf("actual value of constBkg ---------------------------------------->> | %f \n", GetBkgInvMassConst());
+        }else{
+               printf("actual value of m^{0} ------------------------------------------->> | %f \n", GetBkgInvMassNorm());
+               printf("actual value of m^{1} ------------------------------------------->> | %f \n", GetBkgInvMassMean());
+               printf("actual value of m^{2} ------------------------------------------->> | %f \n", GetBkgInvMassSlope());
+               printf("actual value of m^{3} ------------------------------------------->> | %f \n", GetBkgInvMassConst());
+               printf("actual value of m^{4} ------------------------------------------->> | %f \n", GetPolyn4());
+               printf("actual value of m^{5} ------------------------------------------->> | %f \n", GetPolyn5());
+        }
+
        // resolution func (First-First)
        printf("actual value of norm1Gauss (FF) --------------------------------->> | %f \n", GetNormGaus1ResFunc(2));
        printf("actual value of norm2Gauss (FF) --------------------------------->> | %f \n", GetNormGaus2ResFunc(2));
@@ -266,36 +321,37 @@ void AliDielectronBtoJPSItoEleCDFfitFCN::SetResolutionConstants(const Double_t*
 }
 
 //________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m, Int_t type) const 
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m, Double_t pt, Int_t type) const 
 {
         // evaluate likelihood function
-       return fParameters[8]*EvaluateCDFfuncSignalPart(x,m,type) + (1. - fParameters[8])*EvaluateCDFfuncBkgPart(x,m,type);
+       //printf("CDF func ---> x = %f m = %f pt = %f type = %d \n",x,m,pt,type);
+       return fParameters[8]*EvaluateCDFfuncSignalPart(x,m, pt, type) + (1. - fParameters[8])*EvaluateCDFfuncBkgPart(x,m, pt, type);
 }
 
 //_________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m, Int_t type) const
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m, Double_t pt, Int_t type) const
 {
         // evaluate likelihood function
-       return EvaluateCDFfunc(x,m,type);
+       return EvaluateCDFfunc(x,m,pt, type);
 }
 
 //_________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m, Int_t type) const 
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m, Double_t pt, Int_t type) const 
 {
   // evaluate psproper signal  
-  return EvaluateCDFDecayTimeSigDistr(x,type)*(EvaluateCDFInvMassSigDistr(m)/fintmMassSig); 
+  return EvaluateCDFDecayTimeSigDistr(x,pt, type)*(EvaluateCDFInvMassSigDistr(m)/fintmMassSig); 
 }
 
 //_________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x, Int_t type) const
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x, Double_t pt, Int_t type) const
 {
        //
        // Implementation of the Background part of the Likelyhood function
        // 
        Double_t retvalue = 0.;
-       Double_t funBnorm = FunB(x,type);
-       Double_t funPnorm = ResolutionFunc(x,type);
-       retvalue = fParameters[7]*funBnorm + (1. - fParameters[7])*funPnorm;
+       Double_t funBnorm =  fMultivariate ? FunBsaved(x, pt, type) : FunB(x,pt, type)  ;
+        Double_t funPnorm = ResolutionFunc(x, pt, type);
+        retvalue = fParameters[7]*funBnorm + (1. - fParameters[7])*funPnorm;
        return retvalue;
 }
 
@@ -307,9 +363,9 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t
        // It can be either Crystal Ball function or sum of two Landau
        //
        Double_t fitval = 0.;
-
+        // change inv Mass RMS fChangeMass 
        if(fCrystalBallParam){
-               Double_t t = (m-fParameters[9])/fParameters[11]; ;
+               Double_t t = ((m-fParameters[9])/fChangeMass)/fParameters[11]; ;
                if (fParameters[12] < 0) t = -t;
                Double_t absAlpha = TMath::Abs((Double_t)fParameters[12]);
 
@@ -331,14 +387,14 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t
        }
 }
 //_________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunB(Double_t x, Int_t type) const  
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunB(Double_t x, Double_t pt, Int_t type) const  
 {
        //  
        // Parameterisation of the fit function for the x part of the Background
        //
-       Double_t np = 1000.0;
+       Double_t np = 1000.0;
        Double_t sc = 10.;
-       Double_t sigma3 = 1000.; // valore usato nella macro
+       Double_t sigma3 = 1000.;
        Double_t xprime;
        Double_t sum = 0.0;
        Double_t xlow,xupp;
@@ -355,19 +411,20 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunB(Double_t x, Int_t type) const
                xprime = xlow + (i-.5) * step;
                csiMCxprime = CsiMC(xprime);
                xdiff = xprime - x;
-               resolutionxdiff = ResolutionFunc(xdiff, type); // normalized value
+               resolutionxdiff = ResolutionFunc(xdiff, pt, type); // normalized value
                sum += csiMCxprime * resolutionxdiff;
        }
-
+     
        return step * sum ;
 }
+
 //_________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunP(Double_t x, Int_t type) const 
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunP(Double_t x, Double_t pt, Int_t type) const 
 {
        //
        //  Parameterisation of the Prompt part for the x distribution
        //
-       return ResolutionFunc(x,type);
+       return ResolutionFunc(x, pt, type);
 }
 
 
@@ -379,30 +436,36 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::CsiMC(Double_t x) const
        //  for the J/psi coming from Beauty hadrons
        //
        Double_t returnvalue = 0.;
+       
+       if((fhCsiMC->FindBin(x-fShiftTemplate) > 0) && (fhCsiMC->FindBin(x-fShiftTemplate) < fhCsiMC->GetNbinsX()+1))  
+       returnvalue = fhCsiMC->Interpolate(x-fShiftTemplate);
 
-       if((fhCsiMC->FindBin(x) > 0) && (fhCsiMC->FindBin(x) < fhCsiMC->GetNbinsX()+1))  
-       returnvalue = fhCsiMC->Interpolate(x);
 
        return returnvalue;
 }
 
 //_________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m, Int_t type) const 
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m, Double_t pt, Int_t type) const 
 {
        //
        // Return the part of the likelihood function for the background hypothesis
        //
-       return EvaluateCDFDecayTimeBkgDistr(x,type)*(EvaluateCDFInvMassBkgDistr(m)/fintmMassBkg);
-}  
+         Double_t bkgValx = fMultivariate ? EvaluateCDFDecayTimeBkgDistrSaved(x,type,m,pt) : EvaluateCDFDecayTimeBkgDistr(x,type,m,pt);
+        return bkgValx*(EvaluateCDFInvMassBkgDistr(m)/fintmMassBkg);
+}
+  
 
 //_________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x, Int_t type) const 
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x, Int_t type, Double_t m, Double_t pt) const
 {
-       //
-       // it returns the value of the probability to have a given x for the background 
-       //
-
-       Double_t ret = fParameters[0]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*ResolutionFunc(x,type) + fParameters[1]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*FunBkgPos(x,type) +  fParameters[2]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*FunBkgNeg(x,type) + fParameters[3]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*FunBkgSym(x,type);
+        //
+        // it returns the value of the probability to have a given x for the background 
+        // in the pt, m , type correspondent range 
+        //
+       Double_t ret = 0.;
+       if(fMultivariate) 
+       ret = EvaluateCDFDecayTimeBkgDistrDifferential(x,type,m,pt);
+       else ret = fParameters[0]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*ResolutionFunc(x, pt, type) + fParameters[1]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgPos(x, pt,type) +  fParameters[2]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgNeg(x,pt,type) + fParameters[3]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgSym(x, pt,type) + fParameters[46]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgSym1(x,pt,type);
        return ret;
 }
 
@@ -413,18 +476,22 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistr(Double_t
        // it returns the value of the probability to have a given mass for the background
        //
        Double_t value = 0.;
-       value = fParameters[14]*TMath::Exp(-1*(m-fParameters[15])/fParameters[16]) + fParameters[17];  
-       return value;
+       if(fExponentialParam) 
+        value = fParameters[14]*TMath::Exp(-1*(m-fParameters[15])/fParameters[16]) + fParameters[17];  
+        else value = fParameters[14] + fParameters[15]*m + fParameters[16]*m*m + fParameters[17]*m*m*m + fParameters[47]*m*m*m*m + fParameters[48]*m*m*m*m*m; 
+        return value;
 }
 //_________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x,Int_t type) const 
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x, Double_t pt, Int_t type) const 
 {
-       //
+        //
        // exponential with positive slopes for the background part (x)
        //
-       Double_t np = 1000.0;
-       Double_t sc = 10.;      
-       Double_t sigma3 = 1000.; // valore usato nella macro 
+       Double_t np, sc, sigma3;
+       sc = 10.;
+       if(fMultivariate){ np = 10000.0; sigma3 = 5000.;}
+       else{ np = 1000.0; sigma3 = 1000.;}
+
        Double_t xprime;
        Double_t sum = 0.0;
        Double_t xlow,xupp;
@@ -436,22 +503,24 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x,Int_t type) co
 
        for(i=1.0; i<=np/2; i++) {
                xprime = xlow + (i-.5) * step;
-               if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,type));}
+               if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,pt,type));}
                xprime = xupp - (i-.5) * step;
-               if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,type));}
-       }
+               if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,pt,type));}
+               }
 
        return step * sum ;
 }
 //_________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x, Int_t type) const 
-{
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x, Double_t pt, Int_t type) const 
+{ 
        //
        // exponential with negative slopes for the background part (x)
        //
-       Double_t np = 1000.0;
-       Double_t sc = 10.;      
-       Double_t sigma3 = 1000.; 
+       Double_t np, sc, sigma3;
+        sc = 10.;
+        if(fMultivariate){ np = 10000.0;  sigma3 = 5000.;}
+        else{ np = 1000.0; sigma3 = 1000.;}
+       
        Double_t xprime;
        Double_t sum = 0.0;
        Double_t xlow,xupp;
@@ -464,23 +533,25 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x, Int_t type) c
        for(i=1.0; i<=np/2; i++) {
 
                xprime = xlow + (i-.5) * step;
-               if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,type));}
+               if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,pt,type));}
 
                xprime = xupp - (i-.5) * step;
-               if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,type));}
+               if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,pt,type));}
        }
 
        return step * sum ;
 }
 //_________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x, Int_t type) const 
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x, Double_t pt, Int_t type) const 
 {
        //
        // exponential with both positive and negative slopes for the background part (x)
        //
-       Double_t np = 1000.0;
-       Double_t sc = 10.;      
-       Double_t sigma3 = 1000.; 
+       Double_t np, sc, sigma3;
+        sc = 10.;
+        if(fMultivariate){ np = 10000.0; sigma3 = 5000.;}
+        else{ np = 1000.0; sigma3 = 1000.;}
+
        Double_t xprime;
        Double_t sum1 = 0.0;
        Double_t sum2 = 0.0;
@@ -494,22 +565,60 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x, Int_t type) c
        for(i=1.0; i<=np/2; i++) {
 
                xprime = xlow + (i-.5) * step;
-               if (xprime > 0) {sum1 += 0.5 * fParameters[6]*TMath::Exp(-1*xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));}
-               if (xprime < 0) {sum2 += 0.5 * fParameters[6]*TMath::Exp(xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));}
+               if (xprime > 0) {sum1 += 0.5 * fParameters[6]*TMath::Exp(-1*xprime*fParameters[6]) * (ResolutionFunc(xprime-x,pt,type));}
+               if (xprime < 0) {sum2 += 0.5 * fParameters[6]*TMath::Exp(xprime*fParameters[6]) * (ResolutionFunc(xprime-x,pt,type));}
 
                xprime = xupp - (i-.5) * step;
-               if (xprime > 0) {sum1 += 0.5 * fParameters[6]*TMath::Exp(-1*xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));} 
-               if (xprime < 0) {sum2 += 0.5 * fParameters[6]*TMath::Exp(xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));}
+               if (xprime > 0) {sum1 += 0.5 * fParameters[6]*TMath::Exp(-1*xprime*fParameters[6]) * (ResolutionFunc(xprime-x,pt,type));} 
+               if (xprime < 0) {sum2 += 0.5 * fParameters[6]*TMath::Exp(xprime*fParameters[6]) * (ResolutionFunc(xprime-x,pt,type));}
        }
 
        return step*(sum1 + sum2) ;
 }
+
+//_________________________________________________________________________________________________
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgSym1(Double_t x, Double_t pt, Int_t type) const
+{
+       //
+        // exponential with both positive and negative slopes for the background part (x)
+        //
+        Double_t np, sc, sigma3;
+        sc = 10.;
+        if(fMultivariate){ np = 10000.0; sigma3 = 5000.;}
+        else{ np = 1000.0; sigma3 = 1000.;}
+
+       Double_t xprime;
+        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++) {
+
+                xprime = xlow + (i-.5) * step;
+                if (xprime > 0) {sum1 += 0.5 * fParameters[45]*TMath::Exp(-1*xprime*fParameters[45]) * (ResolutionFunc(xprime-x,pt,type));}
+                if (xprime < 0) {sum2 += 0.5 * fParameters[45]*TMath::Exp(xprime*fParameters[45]) * (ResolutionFunc(xprime-x,pt,type));}
+
+                xprime = xupp - (i-.5) * step;
+                if (xprime > 0) {sum1 += 0.5 * fParameters[45]*TMath::Exp(-1*xprime*fParameters[45]) * (ResolutionFunc(xprime-x,pt,type));} 
+                if (xprime < 0) {sum2 += 0.5 * fParameters[45]*TMath::Exp(xprime*fParameters[45]) * (ResolutionFunc(xprime-x,pt,type));}
+        }
+
+        return step*(sum1 + sum2) ;
+}
+
+
 //_________________________________________________________________________________________________
-Double_t AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x, Int_t type) const  
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x, Double_t pt, Int_t type) const  
 {
        //
        // parametrization with 2 gaus
        //
+        x = x/fChangeResolution;
         Int_t index = (2-type)*9;
         Double_t mean1 = fParameters[20+index];
         Double_t mean2 = fParameters[22+index];
@@ -522,15 +631,32 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x, Int_t ty
         Double_t lambda = fParameters[25+index];
         Double_t norm3 = fParameters[26+index];
  
+        if(fMultivariate) // set parameters from matrix
+        {
+          //pt dependence
+          Int_t binPt = -1.;
+          for(int j=0; j<fPtWindows->GetSize()-1; j++) {if(fPtWindows->At(j)<pt && pt<fPtWindows->At(j+1)) binPt = j;}
+          norm1 = fResParams[binPt][type][0];
+          mean1 = fResParams[binPt][type][1];
+          sigma1 = fResParams[binPt][type][2];
+          norm2 = fResParams[binPt][type][3];
+          mean2 = fResParams[binPt][type][4];
+          sigma2 = fResParams[binPt][type][5];
+          alfa = fResParams[binPt][type][6];
+          lambda = fResParams[binPt][type][7];
+          norm3 = fResParams[binPt][type][8];
+        }
+
         Double_t ret = 0.; Double_t fitval = 0; 
         if(TMath::Abs(x)<=alfa) fitval = (lambda-1)/(2*alfa*lambda);
         else  fitval = ((lambda-1)/(2*alfa*lambda))*TMath::Power(alfa,lambda)*(TMath::Power(TMath::Abs(x),-1*lambda));
 
         ret = (norm1/(norm1+norm2+norm3))*((1/(sigma1*TMath::Sqrt(2*TMath::Pi())))*TMath::Exp(-0.5*((x-mean1)/sigma1)*((x-mean1)/sigma1))) + (norm2/(norm1+norm2+norm3))*((1/(sigma2*TMath::Sqrt(2*TMath::Pi())))*TMath::Exp(-0.5*((x-mean2)/sigma2)*((x-mean2)/sigma2))) + (norm3/(norm1+norm2+norm3))*fitval;
 
-        return ret;
+        return ret/fChangeResolution;
 }     
 
+
 //_________________________________________________________________________________________________
 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetCsiMC(Double_t xmin, Double_t xmax, Double_t normalization) 
 {
@@ -542,11 +668,12 @@ TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetCsiMC(Double_t xmin, Double_t xmax,
 }
 
 //__________________________________________________________________________________________________
-TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetResolutionFunc(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
+TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetResolutionFunc(Double_t xmin, Double_t xmax, Double_t normalization, Double_t pt, Int_t type){
        // return the pointer to the resolution function
-       TF1* resFunc = new TF1(Form("resolutionFunc_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFuncf,xmin,xmax,2);
+       TF1* resFunc = new TF1(Form("resolutionFunc_%d",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFuncf,xmin,xmax,3);
         resFunc->SetParameter(0,normalization);
-        resFunc->SetParameter(1,type);
+        resFunc->SetParameter(1,pt);
+        resFunc->SetParameter(2,(Double_t)type);
        resFunc->SetNpx(5000);
         return (TF1*)resFunc->Clone();
 }
@@ -561,22 +688,24 @@ TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetResolutionFuncAllTypes(Double_t xmin
 }
 
 //___________________________________________________________________________________________________
-TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeBkgDistr(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
+TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeBkgDistr(Double_t xmin, Double_t xmax, Double_t normalization, Int_t type, Double_t mass, Double_t pt, Int_t npx){
        // return the pointer to the background x distribution function
-       TF1 *backFunc = new TF1(Form("backFunc_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrFunc,xmin,xmax,2);
+       TF1 *backFunc = new TF1(Form("backFunc_%d_%1.2f_%1.2f",type,mass,pt),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrFunc,xmin,xmax,4);
         backFunc->SetParameter(0,normalization);
-        backFunc->SetParameter(1,type);
-       backFunc->SetNpx(5000);
+        backFunc->SetParameter(1,(Double_t)type);
+        backFunc->SetParameter(2,mass);
+        backFunc->SetParameter(3,pt);
+       backFunc->SetNpx(npx);
         return (TF1*)backFunc->Clone();
 }
 
 //___________________________________________________________________________________________________
 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeBkgDistrAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
         // return the pointer to the background x distribution function
-        TF1 *backFunc = new TF1(Form("backFunc"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrFuncAllTypes,xmin,xmax,1);
-        backFunc->SetParameter(0,normalization);
-        backFunc->SetNpx(5000);
-        return (TF1*)backFunc->Clone();
+        TF1 *backFuncNew = new TF1(Form("backFunc_%f",normalization),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrFuncAllTypes,xmin,xmax,1);
+        backFuncNew->SetParameter(0,normalization);
+        backFuncNew->SetNpx(5000);
+        return (TF1*)backFuncNew->Clone();
 }
 
 //__________________________________________________________________________________________________
@@ -628,11 +757,12 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassTotalDistr(const
 }
 
 //____________________________________________________________________________________________________
-TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeTotalDistr(Double_t xMin, Double_t xMax, Double_t normalization,Double_t type){
+TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeTotalDistr(Double_t xMin, Double_t xMax, Double_t normalization,Double_t pt, Int_t type){
  // return the pointer to the pseudoproper distribution for the background
- TF1 *decayTimeTot = new TF1(Form("decayTimeTot_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistr,xMin,xMax,2);
+ TF1 *decayTimeTot = new TF1(Form("decayTimeTot_%d",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistr,xMin,xMax,3);
  decayTimeTot->SetParameter(0,normalization);
- decayTimeTot->SetParameter(1,type);
+ decayTimeTot->SetParameter(1,pt);
+ decayTimeTot->SetParameter(2,(Double_t)type);
  decayTimeTot->SetNpx(5000);
  return (TF1*)decayTimeTot->Clone();
 }
@@ -643,7 +773,7 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistr(cons
  // evaluate the total pseudoproper distribution for a given candidate's type. par[1] should be the candidate's type.
  Double_t value = 0;
  Double_t xx = x[0];
- value = (fParameters[8]*EvaluateCDFDecayTimeSigDistr(xx,(Int_t)par[1]) + (1-fParameters[8])*EvaluateCDFDecayTimeBkgDistr(xx,(Int_t)par[1]))*par[0];
+ value = (fParameters[8]*EvaluateCDFDecayTimeSigDistr(xx,par[1],(Int_t)par[2]) + (1-fParameters[8])*EvaluateCDFDecayTimeBkgDistr(xx,(Int_t)par[2],par[1]))*par[0];
  return value;
 }
 
@@ -654,7 +784,7 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistrAllTy
  Double_t value = 0;
  Double_t xx = x[0];
 
- value = (fParameters[8]*(fWeightType[2]*EvaluateCDFDecayTimeSigDistr(xx,2)+fWeightType[1]*EvaluateCDFDecayTimeSigDistr(xx,1)+fWeightType[0]*EvaluateCDFDecayTimeSigDistr(xx,0)))+((1-fParameters[8])*(fWeightType[2]*EvaluateCDFDecayTimeBkgDistr(xx,2) + fWeightType[1]*EvaluateCDFDecayTimeBkgDistr(xx,1)+fWeightType[0]*EvaluateCDFDecayTimeBkgDistr(xx,0))); 
+ value = (fParameters[8]*(fWeightType[2]*EvaluateCDFDecayTimeSigDistr(xx,200.,2)+fWeightType[1]*EvaluateCDFDecayTimeSigDistr(xx,200.,1)+fWeightType[0]*EvaluateCDFDecayTimeSigDistr(xx,200.,0)))+((1-fParameters[8])*(fWeightType[2]*EvaluateCDFDecayTimeBkgDistr(xx,2,3.09,200.) + fWeightType[1]*EvaluateCDFDecayTimeBkgDistr(xx,1,3.09,200.)+fWeightType[0]*EvaluateCDFDecayTimeBkgDistr(xx,0,3.09,200.))); 
 
  return value*par[0];
 }
@@ -670,22 +800,175 @@ TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeTotalDistrAllTyp
 
 
 //____________________________________________________________________________________________________
-TF1 * AliDielectronBtoJPSItoEleCDFfitFCN::GetFunB(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
+TF1 * AliDielectronBtoJPSItoEleCDFfitFCN::GetFunB(Double_t xmin, Double_t xmax, Double_t normalization, Double_t pt, Int_t type, Int_t npx){
  // return the pointer to the function that describe secondary jpsi pseudoproper distribution for a given candidate's type
- TF1* funb = new TF1(Form("secondaryJpsiConvolution_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::FunBfunc,xmin,xmax,2);
+ TF1* funb = new TF1(Form("secondaryJpsiConvolution_%d_%1.2f",type,pt),this,&AliDielectronBtoJPSItoEleCDFfitFCN::FunBfunc,xmin,xmax,3);
  funb->SetParameter(0,normalization);
- funb->SetParameter(1,type);
- funb->SetNpx(5000);
+ funb->SetParameter(1,pt);
+ funb->SetParameter(2,(Double_t)type);
+ funb->SetNpx(npx);
  return (TF1*)funb->Clone();
  }
 
 //____________________________________________________________________________________________________
 TF1 * AliDielectronBtoJPSItoEleCDFfitFCN::GetFunBAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
-// return the pointer to the function that describe secondary jpsi pseudoproper distribution for all types
+ // return the pointer to the function that describe secondary jpsi pseudoproper distribution for all types
  TF1* funb = new TF1(Form("secondaryJpsiConvolution"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::FunBfuncAllTypes,xmin,xmax,1);
  funb->SetParameter(0,normalization);
  funb->SetNpx(5000);
  return (TF1*)funb->Clone();
  }
 
+//
+// methods below are needed to perform the multivariate fit (pt,mass,type); this can be enabled 
+// by the boolean fMultivariate; if functions to describe pseudoproper
+// decay lenght in pt,m, type bins have been
+// already computed, they can be loaded from the file 
+// switching-on the option fLoadFunction (this to avoid the
+// calculation of convolutions and speed-up the likelihood fit)
+//
+
+//____________________________________________________________________________________________________
+void AliDielectronBtoJPSItoEleCDFfitFCN::SetBackgroundSpecificParameters(Int_t pt, Int_t mb, Int_t tp){
+  //
+  // methods to set specific background parameters in bins(pt, inv. mass, type)
+  //
+  for(int j=0; j<4;j++) fParameters[j]=fBkgParams[pt][mb][tp][j];
+  for(int k=5; k<8;k++) fParameters[k-1]=fBkgParams[pt][mb][tp][k];
+  fParameters[45] = fBkgParams[pt][mb][tp][8];
+  fParameters[46] = fBkgParams[pt][mb][tp][4];
+  return;
+}
+
+//_______________________________________________________________________________________________________
+void AliDielectronBtoJPSItoEleCDFfitFCN::InitializeFunctions(Int_t ptSize, Int_t massSize){
+  //
+  // initialize pointers to save functions for the  multivariate fit
+  //
+  fFunBSaved = new TF1**[ptSize];
+  for(int kpt=0; kpt<ptSize; kpt++) fFunBSaved[kpt] = new TF1*[3]; // type
+  fFunBkgSaved = new TF1***[ptSize];
+  for(int kpt=0; kpt<ptSize; kpt++){ fFunBkgSaved[kpt] = new TF1**[massSize];
+  for(int ks = 0; ks<massSize; ks++) fFunBkgSaved[kpt][ks] = new TF1*[3];
+  for(int kk=0; kk<3;kk++) {
+     fFunBSaved[kpt][kk] = new TF1();
+     for(int kk1=0; kk1<massSize;kk1++){ 
+     fFunBkgSaved[kpt][kk1][kk] = new TF1();
+     }
+    }
+  }
+
+  // to extrapolate the function under the signal region
+  fWeights = new Double_t**[massSize - 1]; // mass
+  for(int km =0; km < (massSize - 1); km++) {fWeights[km] = new Double_t*[ptSize];
+        for(int kpt =0; kpt<ptSize; kpt++) fWeights[km][kpt] = new Double_t[3];
+  } // pt
+ return;
+}
+
+//_________________________________________________________________________________________________
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBsaved(Double_t x, Double_t pt, Int_t type) const
+{
+        //
+        //   functions to describe non-prompt J/psi x distributions
+        //
+        Double_t returnvalue = 0.;
+        Int_t binPt = -1;
+        for(int j=0; j<fPtWindows->GetSize()-1; j++) {if(fPtWindows->At(j)<pt && pt<fPtWindows->At(j+1)) binPt = j;}
+        returnvalue = fFunBSaved[binPt][type]->Eval(x);
+        return returnvalue;
+}
+
+//________________________________________________________________________________________________
+void AliDielectronBtoJPSItoEleCDFfitFCN::SetFunctionsSaved(Int_t npxFunB, Int_t npxFunBkg, Double_t funBLimits, Double_t funBkgLimits, Int_t signalRegion){
+       //
+       // save functions for the multivariate fit 
+       //
+       if(!fMultivariate)
+       {AliInfo("Warning: fMultivariate is kFALSE! Functions are not saved! \n"); return;}
+       SetExtrapolationRegion(signalRegion);
+
+                     for(int tp=0;tp<3;tp++)  // type
+                      {
+                        // pt
+                        for(int pt=0; pt<fPtWindows->GetSize()-1;pt++){
+                        if(fResParams) SetResolutionConstants(fResParams[pt][tp],tp);
+                        SetFunBFunction(tp,pt,GetFunB(-1.*funBLimits,funBLimits,1.,(fPtWindows->At(pt) + (fPtWindows->At(pt+1)-fPtWindows->At(pt))/2.),tp,npxFunB));
+                        }
+                      }
+
+       AliInfo("+++++++  Pseudoproper-decay-length function for secondary J/psi saved  ++++++ \n");
+
+       if(!fLoadFunctions){
+       for(int ij = 0; ij<fMassWindows->GetSize()-1;ij++){
+       if(ij == signalRegion) continue;
+
+       Int_t mbin = (ij > signalRegion) ?  ij-1 : ij;
+       for(int tp=0;tp<3;tp++)  {
+       for(int pt =0; pt<fPtWindows->GetSize()-1; pt++){
+         if(fBkgParams) SetBackgroundSpecificParameters(pt,mbin,tp);
+         SetBkgFunction(ij, tp, pt, GetEvaluateCDFDecayTimeBkgDistr(-1.*funBkgLimits,funBkgLimits,1.,tp,(fMassWindows->At(ij) + (fMassWindows->At(ij+1)-fMassWindows->At(ij))/2.),(fPtWindows->At(pt) + (fPtWindows->At(pt+1)-fPtWindows->At(pt))/2.),npxFunBkg));
+
+               }
+
+       }
+
+       }
+       AliInfo("+++++++  Pseudoproper-decay-length function for background saved  +++++++++++ \n");
+
+       } // loadFunctions
+       // evaluate under signal
+       for(int tp=0;tp<3;tp++)
+                       {
+               for(int pt =0; pt<fPtWindows->GetSize()-1; pt++){
+               SetBkgFunction(signalRegion, tp, pt, GetEvaluateCDFDecayTimeBkgDistr(-1.*funBkgLimits,funBkgLimits,1.,tp,(fMassWindows->At(signalRegion) + (fMassWindows->At(signalRegion+1)-fMassWindows->At(signalRegion))/2.),(fPtWindows->At(pt) + (fPtWindows->At(pt+1)-fPtWindows->At(pt))/2.),npxFunBkg));
+               }
+   
+       }
+       // save functions
+       TFile func("functions.root","RECREATE");
+       for(int kpt =0; kpt<fPtWindows->GetSize()-1; kpt++){
+       for(int ss=0; ss<3;ss++) {fFunBSaved[kpt][ss]->Write();
+       for(int kk=0; kk<fMassWindows->GetSize()-1; kk++) fFunBkgSaved[kpt][kk][ss]->Write();}}
+       return;
+}
+
+//_________________________________________________________________________________________________
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrDifferential(Double_t x, Int_t type, Double_t m, Double_t pt) const
+{
+        //
+        // it returns the value of the probability to have a given x for the background 
+        // in the pt, m , type correspondent range 
+        //
+        Int_t binPt = -1;
+        for(int j=0; j<fPtWindows->GetSize()-1; j++)
+        {if(fPtWindows->At(j)<pt && pt<fPtWindows->At(j+1)) binPt = j;}
+        Bool_t isSignal = (fMassWindows->At(fSignalBinForExtrapolation)<m && m<fMassWindows->At(fSignalBinForExtrapolation+1));
+        Double_t ret = 0.;
+        if(!isSignal)
+        ret = fParameters[0]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*ResolutionFunc(x, pt, type) + fParameters[1]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgPos(x, pt,type) +  fParameters[2]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgNeg(x,pt,type) + fParameters[3]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgSym(x, pt,type) + fParameters[46]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgSym1(x,pt,type);
+       else{
+       for(int k=0; k<fMassWindows->GetSize()-2;k++) {
+        Int_t mbin = (k > (fSignalBinForExtrapolation-1)) ?  k+1 : k;
+        ret +=  fWeights[k][binPt][type]*EvaluateCDFDecayTimeBkgDistrSaved(x,type,(fMassWindows->At(mbin) + (fMassWindows->At(mbin+1)-fMassWindows->At(mbin))/2.),pt);}
+
+       }
+        return ret;
+}
+
+//_________________________________________________________________________________________________
+Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrSaved(Double_t x, Int_t type, Double_t m, Double_t pt) const
+{
+        //
+        // it returns the value of the probability to have a given x for the background 
+        // in the pt, m , type correspondent range 
+        //
+        Double_t returnvalue = 0.;
+        Int_t binM = -1.;
+        for(int j=0; j<fMassWindows->GetSize()-1; j++) {if(fMassWindows->At(j)<m && m<fMassWindows->At(j+1)) binM = j;}
+        Int_t binPt = -1;
+        for(int j=0; j<fPtWindows->GetSize()-1; j++) {if(fPtWindows->At(j)<pt && pt<fPtWindows->At(j+1)) binPt = j;}
+        returnvalue = fFunBkgSaved[binPt][binM][type]->Eval(x);
+        return returnvalue;
+}
 
index 4054f68..1a7aceb 100644 (file)
 \r
 #include <TNamed.h>\r
 #include <TDatabasePDG.h>\r
-#include "TH1F.h"\r
+#include <TH1F.h>\r
+\r
 \r
 class TRandom3;\r
 class TF1;\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 AliDielectronBtoJPSItoEleCDFfitFCN : public TNamed {\r
        public:\r
@@ -39,15 +31,18 @@ class AliDielectronBtoJPSItoEleCDFfitFCN : public TNamed {
                virtual ~AliDielectronBtoJPSItoEleCDFfitFCN();\r
 \r
                Double_t EvaluateLikelihood(const Double_t* pseudoproperdecaytime,\r
-                               const Double_t* invariantmass, const Int_t* type, const Int_t ncand) const;\r
+                               const Double_t* invariantmass, const Double_t* pt, const Int_t* type, const Int_t ncand) const;\r
 \r
+                // getters for parameters\r
                Double_t GetResWeight()                    const { return fParameters[0]; }\r
                Double_t GetFPlus()                        const { return fParameters[1]; }\r
                Double_t GetFMinus()                       const { return fParameters[2]; }\r
                Double_t GetFSym()                         const { return fParameters[3]; } \r
+               Double_t GetFSym1()                        const { return fParameters[46]; }\r
                Double_t GetLamPlus()                      const { return fParameters[4]; }\r
                Double_t GetLamMinus()                     const { return fParameters[5]; }\r
                Double_t GetLamSym()                       const { return fParameters[6]; }\r
+               Double_t GetLamSym1()                      const { return fParameters[45]; }\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
@@ -67,39 +62,40 @@ class AliDielectronBtoJPSItoEleCDFfitFCN : public TNamed {
                 Double_t GetResSigma1(Int_t type)          const { return fParameters[21+(2-type)*9]; }\r
                 Double_t GetResMean2(Int_t type)           const { return fParameters[22+(2-type)*9]; }\r
                 Double_t GetResSigma2(Int_t type)          const { return fParameters[23+(2-type)*9]; }\r
-                \r
                 Double_t GetResAlfa(Int_t type)            const { return fParameters[24+(2-type)*9]; }   \r
                 Double_t GetResLambda(Int_t type)          const { return fParameters[25+(2-type)*9]; } \r
                 Double_t GetResNormExp(Int_t type)         const { return fParameters[26+(2-type)*9]; }\r
-\r
+                Double_t GetPolyn4()                       const { return fParameters[47]; }\r
+                Double_t GetPolyn5()                       const { return fParameters[48]; }\r
                 Bool_t GetCrystalBallParam()               const { return fCrystalBallParam; }\r
+                Bool_t GetExponentialParam()               const { return fExponentialParam; }\r
                TH1F * GetCsiMcHisto()                     const { return fhCsiMC; }\r
                 Double_t GetResWeight(Int_t iW)            const { return fWeightType[iW]; }\r
-               \r
-               Double_t* GetParameters()                         { return fParameters;}\r
 \r
                // return pointer to likelihood functions  \r
                TF1* GetCsiMC(Double_t xmin, Double_t xmax,Double_t normalization);\r
-               TF1* GetResolutionFunc(Double_t xmin, Double_t xmax,Double_t normalization, Double_t type=2);\r
+               TF1* GetResolutionFunc(Double_t xmin, Double_t xmax,Double_t normalization, Double_t pt, Int_t type=2);\r
                TF1* GetResolutionFuncAllTypes(Double_t xmin, Double_t xmax,Double_t normalization);\r
-               TF1* GetFunB(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type=2);\r
+               TF1* GetFunB(Double_t xmin, Double_t xmax, Double_t normalization, Double_t pt, Int_t type=2, Int_t npx = 5000);\r
                 TF1* GetFunBAllTypes(Double_t xmin, Double_t xmax, Double_t normalization);\r
-                TF1* GetEvaluateCDFDecayTimeBkgDistr(Double_t xmin, Double_t xmax, Double_t normalization,Double_t type = 2);\r
+                TF1* GetEvaluateCDFDecayTimeBkgDistr(Double_t xmin, Double_t xmax, Double_t normalization, Int_t type = 2, Double_t mass = 3.09, Double_t pt = 200.,Int_t npx = 5000);\r
                TF1* GetEvaluateCDFDecayTimeBkgDistrAllTypes(Double_t xmin, Double_t xmax, Double_t normalization);\r
                 TF1* GetEvaluateCDFDecayTimeSigDistr(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type);\r
                 TF1* GetEvaluateCDFInvMassBkgDistr(Double_t mMin, Double_t mMax, Double_t normalization);\r
                 TF1* GetEvaluateCDFInvMassSigDistr(Double_t mMin, Double_t mMax, Double_t normalization);\r
                 TF1* GetEvaluateCDFInvMassTotalDistr(Double_t mMin, Double_t mMax, Double_t normalization);\r
-                TF1* GetEvaluateCDFDecayTimeTotalDistr(Double_t xMin, Double_t xMax, Double_t normalization, Double_t type=2);\r
+                TF1* GetEvaluateCDFDecayTimeTotalDistr(Double_t xMin, Double_t xMax, Double_t normalization,Double_t pt = 200., Int_t type=2);\r
                 TF1 *GetEvaluateCDFDecayTimeTotalDistrAllTypes(Double_t xMin, Double_t xMax, Double_t normalization);\r
 \r
                void SetResWeight(Double_t resWgt) {fParameters[0] = resWgt;}\r
                void SetFPlus(Double_t plus) {fParameters[1] = plus;}\r
                void SetFMinus(Double_t minus) {fParameters[2]  = minus;}\r
                void SetFSym(Double_t sym) {fParameters[3] = sym;}\r
+               void SetFSym1(Double_t sym) {fParameters[46] = sym;}\r
                void SetLamPlus(Double_t lamplus) {fParameters[4] = lamplus;}\r
                void SetLamMinus(Double_t lamminus) {fParameters[5] = lamminus;}\r
                void SetLamSym(Double_t lamsym) {fParameters[6] = lamsym;}\r
+               void SetLamSym1(Double_t lamsym) {fParameters[45] = lamsym;}\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
@@ -111,32 +107,56 @@ class AliDielectronBtoJPSItoEleCDFfitFCN : public TNamed {
                void SetBkgInvMassMean(Double_t BkgInvMassMean) {fParameters[15] = BkgInvMassMean;}\r
                void SetBkgInvMassSlope(Double_t BkgInvMassSlope) {fParameters[16] = BkgInvMassSlope;}\r
                void SetBkgInvMassConst(Double_t BkgInvMassConst) {fParameters[17] = BkgInvMassConst;}\r
+                void SetBkgInvMassPolyn4(Double_t coeffPol4) {fParameters[47] = coeffPol4;}\r
+                void SetBkgInvMassPolyn5(Double_t coeffPol5) {fParameters[48] = coeffPol5;}\r
                void SetNormGaus1ResFunc(Double_t norm1) {fParameters[18] = norm1;}\r
                void SetNormGaus2ResFunc(Double_t norm2) {fParameters[19] = norm2;}\r
                void SetAllParameters(const Double_t* parameters);\r
                void SetIntegralMassSig(Double_t integral) { fintmMassSig = integral; }\r
                void SetIntegralMassBkg(Double_t integral) { fintmMassBkg = integral; }\r
                void SetCsiMC(const TH1F* MCtemplate) {fhCsiMC = (TH1F*)MCtemplate->Clone("fhCsiMC");}\r
-\r
+                void SetTemplateShift(Double_t shift = 0.){fShiftTemplate = shift;}\r
\r
                void SetResolutionConstants(const Double_t* resolutionConst, Int_t type);\r
                void SetMassWndHigh(Double_t limit) { fMassWndHigh = TDatabasePDG::Instance()->GetParticle(443)->Mass() + limit ;}\r
                void SetMassWndLow(Double_t limit) { fMassWndLow = TDatabasePDG::Instance()->GetParticle(443)->Mass() - limit ;}\r
                void SetCrystalBallFunction(Bool_t okCB) {fCrystalBallParam = okCB;}\r
+               void SetExponentialFunction(Bool_t okExp) {fExponentialParam = okExp;}\r
  \r
                 void SetWeightType(Double_t wFF, Double_t wFS, Double_t wSS) {fWeightType[0]= wSS; fWeightType[1]= wFS; fWeightType[2]= wFF;}\r
-               void ComputeMassIntegral(); \r
 \r
+                void SetChangeResolution(Double_t change){fChangeResolution = change;}   \r
+                void SetChangeMass(Double_t change){fChangeMass = change;}   \r
+               void ComputeMassIntegral(); \r
                void ReadMCtemplates(Int_t BinNum);\r
-\r
                void PrintStatus();\r
 \r
+               //specific methods for multi-variational fit\r
+               void SetBkgWeights(Double_t ***bkgWgt){for(int kpt=0; kpt<(fPtWindows->GetSize()-1); kpt++){for(int k=0; k<(fMassWindows->GetSize()-2); k++) { for(int ktype=0; ktype<3; ktype++) fWeights[k][kpt][ktype]=bkgWgt[k][kpt][ktype]; }}}\r
\r
+               void SetBkgFunction(Int_t massRange,Int_t type, Int_t ptB, TF1 *histBkg){\r
+                fFunBkgSaved[ptB][massRange][type] = histBkg;\r
+                }\r
+                TF1 *GetBkgFunction(Int_t massRange, Int_t ptB, Int_t type) const {return fFunBkgSaved[ptB][massRange][type];} \r
+                void SetFunBFunction(Int_t type, Int_t ptB, TF1 *histSec) { fFunBSaved[ptB][type] = histSec; }\r
+                void SetBackgroundSpecificParameters(Int_t pt, Int_t mb, Int_t tp);\r
+                void SetExtrapolationRegion(Int_t extrRegion){fSignalBinForExtrapolation = extrRegion;}                                           void SetLoadFunction(Bool_t loadFunc) { fLoadFunctions = loadFunc;}\r
+                void SetMultivariateFit(Bool_t multVar) { fMultivariate = multVar;}\r
+                Bool_t GetMultivariate() const { return fMultivariate;} \r
+                void SetFunctionsSaved(Int_t npxFunB=5000, Int_t npxFunBkg=5000, Double_t funBLimits = 20000., Double_t funBkgLimits = 40000., Int_t signalRegion=2);\r
+                void SetResParams(Double_t ***pars){ fResParams = pars;}\r
+                void SetBkgParams(Float_t ****pars){ fBkgParams = pars;}\r
+                void SetMassWindows(TArrayD *msWnd){ fMassWindows = msWnd;}\r
+                void SetPtWindows(TArrayD *ptWnd){ fPtWindows = ptWnd;}\r
+                void InitializeFunctions(Int_t ptSize, Int_t massSize);\r
                Double_t EvaluateCDFInvMassSigDistr(Double_t m) const ;\r
                Double_t EvaluateCDFInvMassBkgDistr(Double_t m) const;\r
-               Double_t ResolutionFunc(Double_t x, Int_t type) const;\r
-               Double_t EvaluateCDFDecayTimeBkgDistr(Double_t x, Int_t type) const ;\r
+               Double_t ResolutionFunc(Double_t x, Double_t pt, Int_t type) const;\r
+               Double_t EvaluateCDFDecayTimeBkgDistr(Double_t x, Int_t type, Double_t m=3.09, Double_t pt=200.) const ;\r
+       \r
 \r
        private:  \r
-               Double_t fParameters[45];        /*  par[0]  = weightRes;                \r
+               Double_t fParameters[49];        /*  par[0]  = weightRes;                \r
                                                     par[1]  = fPos;\r
                                                     par[2]  = fNeg;\r
                                                     par[3]  = fSym\r
@@ -180,35 +200,56 @@ class AliDielectronBtoJPSItoEleCDFfitFCN : public TNamed {
                                                      par[41] = fSigma2ResFunc;\r
                                                      par[42] = fResAlfa; \r
                                                      par[43] = fResLambda;\r
-                                                     par[44] = fResNormExp;\r
-                                                     */\r
-\r
-               Double_t fFPlus;                     // parameters of the log-likelihood function\r
-               Double_t fFMinus;                    // Slopes of the x distributions of the background \r
+                                                     par[44] = fResNormExp;   \r
+                                                    ------------------------------> parameters to describe x and m à la CDF\r
+                                                    // additional parameters below (added for PbPb analysis): if not used should be fixed\r
+                                                    // by the FitHandler pointer\r
+                                                     par[45] = fOneOvLamSym1; //  additional parameter for background\r
+                                                    par[46] = fSym1;         //  additional parameter for background\r
+                                                     par[47] = fPolyn4; //additional parameter for inv. mass background\r
+                                                     par[48] = fPolyn5; //additional parameter for inv. mass background */\r
+                                                  \r
+               Double_t fFPlus;                     // parameter 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 fintmMassSig;               // integral of invariant mass distribution for the signal\r
                Double_t fintmMassBkg;               // integral of invariant mass distribution for the bkg\r
 \r
                TH1F *fhCsiMC;                       // X distribution used as MC template for JPSI from B\r
+\r
+                Double_t fShiftTemplate;             // to shift the MC template\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
                 Double_t fWeightType[3];             // vector with weights of candidates types (used to draw functions)            \r
-                ////\r
-\r
-               Double_t EvaluateCDFfunc(Double_t x, Double_t m, Int_t type) const ;\r
-               Double_t EvaluateCDFfuncNorm(Double_t x, Double_t m, Int_t type) const ;\r
+                Double_t fChangeResolution;         // constant to change the RMS of the resolution function\r
+                Double_t fChangeMass;               // constant to change the RMS of the signal mass function\r
+                \r
+               // data members for multivariate fit\r
+               Double_t ***fWeights;                // weights to interpolate bkg shape in pt, mass, type bins under signal region\r
+                Bool_t fLoadFunctions;              // boolean to load functions saved \r
+                Bool_t fMultivariate;               // switch-on multivariate fit\r
+                TF1 ***fFunBSaved;                  // pointers to save functions for x of non-prompt J/psi in pt, mass, type bins\r
+                TF1 ****fFunBkgSaved;               // pointers to save functions for x of background in pt, mass, type bins\r
+                Double_t ***fResParams;                     // resolution function parameters in pt and type bins\r
+                Float_t ****fBkgParams;              // x background parameters in pt, mass, type bins\r
+                TArrayD *fMassWindows;              // limits for invariant mass bins\r
+                TArrayD *fPtWindows;                // limits for pt bins\r
+               Bool_t fExponentialParam;            // Boolean to switch to Exponential parameterisation\r
+                Double_t fSignalBinForExtrapolation; // inv. mass region in which extrapolate the background shape\r
+                                                    // starting from background shapes in the near inv. mass regions\r
 \r
-               ////\r
 \r
-               Double_t EvaluateCDFfuncSignalPart(Double_t x, Double_t m, Int_t type) const ;      // Signal part \r
-               Double_t EvaluateCDFDecayTimeSigDistr(Double_t x, Int_t type) const ;\r
-               Double_t EvaluateCDFDecayTimeSigDistrFunc(const Double_t* x, const Double_t *par) const { return par[0]*EvaluateCDFDecayTimeSigDistr(x[0],(Int_t)par[1]);}\r
+               Double_t EvaluateCDFfunc(Double_t x, Double_t m, Double_t pt, Int_t type) const ;\r
+               Double_t EvaluateCDFfuncNorm(Double_t x, Double_t m, Double_t pt, Int_t type) const ;\r
+               Double_t EvaluateCDFfuncSignalPart(Double_t x, Double_t m, Double_t pt, Int_t type) const ;      // Signal part \r
+               Double_t EvaluateCDFDecayTimeSigDistr(Double_t x, Double_t pt, Int_t type) const ;\r
+               Double_t EvaluateCDFDecayTimeSigDistrFunc(const Double_t* x, const Double_t *par) const { return par[0]*EvaluateCDFDecayTimeSigDistr(x[0],par[1],(Int_t)par[2]);}\r
                 Double_t EvaluateCDFInvMassSigDistrFunc(const Double_t* x, const Double_t *par) const {return par[0]*EvaluateCDFInvMassSigDistr(x[0])/fintmMassSig;}\r
-               Double_t EvaluateCDFfuncBkgPart(Double_t x,Double_t m,Int_t type) const ;          // Background part\r
-               Double_t EvaluateCDFDecayTimeBkgDistrFunc(const Double_t* x, const Double_t *par) const { return EvaluateCDFDecayTimeBkgDistr(x[0],(Int_t)par[1])*par[0];}\r
+               Double_t EvaluateCDFfuncBkgPart(Double_t x,Double_t m, Double_t pt, Int_t type) const ;          // Background part\r
+               Double_t EvaluateCDFDecayTimeBkgDistrFunc(const Double_t* x, const Double_t *par) const { return EvaluateCDFDecayTimeBkgDistr(x[0],(Int_t)par[1],par[2],par[3])*par[0];}\r
                 Double_t EvaluateCDFDecayTimeBkgDistrFuncAllTypes(const Double_t* x, const Double_t *par) const {return (fWeightType[2]*EvaluateCDFDecayTimeBkgDistr(x[0],2)+fWeightType[1]*EvaluateCDFDecayTimeBkgDistr(x[0],1)+fWeightType[0]*EvaluateCDFDecayTimeBkgDistr(x[0],0))*par[0];}\r
                 Double_t EvaluateCDFInvMassBkgDistrFunc(const Double_t* x, const Double_t *par) const {return par[0]*EvaluateCDFInvMassBkgDistr(x[0])/fintmMassBkg;} \r
                   \r
@@ -217,18 +258,23 @@ class AliDielectronBtoJPSItoEleCDFfitFCN : public TNamed {
                 ////\r
                 Double_t EvaluateCDFDecayTimeTotalDistrAllTypes(const Double_t* x, const Double_t *par) const;\r
 \r
-               Double_t FunB(Double_t x, Int_t type) const;\r
-               Double_t FunBfunc(const Double_t *x, const Double_t *par) const {return FunB(x[0],(Int_t)par[1])*par[0];}\r
-                Double_t FunBfuncAllTypes(const Double_t *x, const Double_t *par) const {return (fWeightType[2]*FunB(x[0],2)+fWeightType[1]*FunB(x[0],1)+fWeightType[0]*FunB(x[0],0))*par[0];}\r
-                Double_t FunP(Double_t x, Int_t type) const ;\r
+               Double_t FunB(Double_t x, Double_t pt, Int_t type) const;\r
+               Double_t FunBfunc(const Double_t *x, const Double_t *par) const {return FunB(x[0],par[1],(Int_t)par[2])*par[0];}\r
+                Double_t FunBfuncAllTypes(const Double_t *x, const Double_t *par) const {return (fWeightType[2]*FunB(x[0],200.,2)+fWeightType[1]*FunB(x[0],200.,1)+fWeightType[0]*FunB(x[0],200.,0))*par[0];}\r
+                Double_t FunP(Double_t x, Double_t pt,Int_t type) const ;\r
                Double_t CsiMC(Double_t x) const;\r
                Double_t CsiMCfunc(const Double_t* x, const Double_t *par) const {  return CsiMC(x[0])*par[0];}\r
-               Double_t FunBkgPos(Double_t x, Int_t type) const ;\r
-               Double_t FunBkgNeg(Double_t x, Int_t type) const ;\r
-               Double_t FunBkgSym(Double_t x, Int_t type) const ;\r
-               Double_t ResolutionFuncf(const Double_t* x, const Double_t *par) const { return ResolutionFunc(x[0],(Int_t)par[1])*par[0];}\r
-                Double_t ResolutionFuncAllTypes(const Double_t* x, const Double_t *par) const { return (fWeightType[2]*ResolutionFunc(x[0],2)+fWeightType[1]*ResolutionFunc(x[0],1)+fWeightType[0]*ResolutionFunc(x[0],0))*par[0]; }                 \r
\r
+               Double_t FunBkgPos(Double_t x, Double_t pt, Int_t type) const ;\r
+               Double_t FunBkgNeg(Double_t x, Double_t pt, Int_t type) const ;\r
+               Double_t FunBkgSym(Double_t x, Double_t pt, Int_t type) const ;\r
+               Double_t FunBkgSym1(Double_t x, Double_t pt, Int_t type) const ;\r
+               Double_t ResolutionFuncf(const Double_t* x, const Double_t *par) const { return ResolutionFunc(x[0],par[1],(Int_t)par[2])*par[0];}\r
+                Double_t ResolutionFuncAllTypes(const Double_t* x, const Double_t *par) const { return (fWeightType[2]*ResolutionFunc(x[0],200.,2)+fWeightType[1]*ResolutionFunc(x[0],200.,1)+fWeightType[0]*ResolutionFunc(x[0],200.,0))*par[0]; }                 \r
+\r
+                // private methods for multivariate fit \r
+                Double_t EvaluateCDFDecayTimeBkgDistrSaved(Double_t x, Int_t type, Double_t m=3.09, Double_t pt = 200.) const ;\r
+               Double_t EvaluateCDFDecayTimeBkgDistrDifferential(Double_t x, Int_t type, Double_t m=3.09, Double_t pt = 200.) const;\r
+                Double_t FunBsaved(Double_t x, Double_t pt, Int_t type) const;\r
 \r
                ClassDef (AliDielectronBtoJPSItoEleCDFfitFCN,1);         // Unbinned log-likelihood fit \r
 \r
index 4079622..15f0d30 100644 (file)
@@ -345,7 +345,7 @@ Double_t AliDielectronBtoJPSItoEleCDFfitFCNfitter::CDFResolutionFunction(const D
 //   resParam[i] = par[i];
 //  }
  fFCN->SetResolutionConstants(par,(Int_t)par[9]);
- return par[10]*fFCN->ResolutionFunc(x[0],(Int_t)(par[9]));
+ return par[10]*fFCN->ResolutionFunc(x[0],-1.,(Int_t)(par[9]));
 }
 //__________________________________________________________________________________________
 Double_t AliDielectronBtoJPSItoEleCDFfitFCNfitter::PsProperBackFitFunc(const Double_t* x, const Double_t* par)
index b8cb278..03adbae 100644 (file)
@@ -46,11 +46,12 @@ void CDFFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t i
 
 //_________________________________________________________________________________________________
 AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler():
-       fIsParamFixed(45),
+       fIsParamFixed(49),
        fPrintStatus(kFALSE),
        fUp(0),
        fX(0x0),
        fM(0x0),
+       fPt(0x0),
        fType(0x0),
         fLikely(0x0),
        fNcand(0),
@@ -62,16 +63,16 @@ AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler()
        //
        // default constructor
        //
-       for (Int_t i=0; i<45; ++i) fParamStartValues[i]=0;
 }
 //_________________________________________________________________________________________________
 AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(Double_t* decaytime, 
-               Double_t* invariantmass, Int_t *type, Int_t ncand) :
-       fIsParamFixed(45),
+               Double_t* invariantmass, Double_t *pt,Int_t *type, Int_t ncand) :
+       fIsParamFixed(49),
        fPrintStatus(kFALSE),
        fUp(0),
        fX(decaytime),
        fM(invariantmass),
+       fPt(pt),
         fType(type),
        fLikely(0x0),
        fNcand(ncand),
@@ -83,7 +84,6 @@ AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(D
        //
        // constructor
        //
-        for (Int_t i=0; i<45; ++i) fParamStartValues[i]=0;
        AliInfo("\n+++\n+++ Minimization object AliDielectronBtoJPSItoEleCDFfitHandler created\n+++\n");
        fLikely = new AliDielectronBtoJPSItoEleCDFfitFCN();
        AliInfo("\n+++\n+++ CDF fit function object AliDielectronBtoJPSItoEleCDFfitFCN created\n+++\n");
@@ -132,6 +132,10 @@ AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(D
         AliInfo("Parameter 42 ----> fAlfaRes (SS)");
         AliInfo("Parameter 43 ----> fLambdaRes (SS)");
         AliInfo("Parameter 44 ----> fNormResExp (SS)");
+        AliInfo("Parameter 45 ----> fOneOvLamSym1 (additional parameter)");
+        AliInfo("Parameter 46 ----> fSym1 (additional parameter)");
+        AliInfo("Parameter 47 ----> fPolyn4 (additional parameter)");
+        AliInfo("Parameter 48 ----> fPolyn5 (additional parameter)");
 
        AliInfo(Form("\n+++\n+++ Number of candidates ---> %d\n+++\n ", ncand));
 }
@@ -147,6 +151,7 @@ AliDielectronBtoJPSItoEleCDFfitHandler& AliDielectronBtoJPSItoEleCDFfitHandler::
                fUp           = c.fUp;
                fX            = c.fX;
                fM            = c.fM;
+               fPt           = c.fPt;
                fType         = c.fType;
                 fLikely       = c.fLikely;
                fNcand        = c.fNcand;
@@ -165,6 +170,7 @@ AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(c
        fUp(c.fUp),
        fX(c.fX),
        fM(c.fM),
+       fPt(c.fPt),
        fType(c.fType),
         fLikely(c.fLikely),
        fNcand(c.fNcand),
@@ -176,7 +182,6 @@ AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(c
        //
        // Copy Constructor
        //
-        for (Int_t i=0; i<45; ++i) fParamStartValues[i]=c.fParamStartValues[i];
 }
 //_______________________________________________________________________________________
 AliDielectronBtoJPSItoEleCDFfitHandler::~AliDielectronBtoJPSItoEleCDFfitHandler()
@@ -192,9 +197,8 @@ Int_t AliDielectronBtoJPSItoEleCDFfitHandler::DoMinimization(Int_t step)
        //
        // performs the minimization
        //
-       if(step == 0){ 
-               //fitter = TVirtualFitter::Fitter(this,20);
-               fitter = (TFitter*)TVirtualFitter::Fitter(this,45);
+       if(step == 0 || !fitter){ 
+               fitter = (TFitter*)TVirtualFitter::Fitter(this,49);
                fitter->SetFCN(CDFFunction);
                fitter->SetParameter(0,"fWeightRes",fParamStartValues[0], 1.e-08, 0., 1.e+06);
                fitter->SetParameter(1,"fPos",fParamStartValues[1], 1.e-08, 0.,1.e+06);
@@ -210,10 +214,10 @@ Int_t AliDielectronBtoJPSItoEleCDFfitHandler::DoMinimization(Int_t step)
                fitter->SetParameter(11,"fSigma",fParamStartValues[11], 1.e-08, 0., 1.e+04);
                fitter->SetParameter(12,"fAlpha",fParamStartValues[12], 1.e-08, 0., 1.e+04);
                fitter->SetParameter(13,"fNorm",fParamStartValues[13], 1.e-08, 0., 1.e+04);
-               fitter->SetParameter(14,"fBkgNorm",fParamStartValues[14], 1.e-08, 0., 1.e+04);
-               fitter->SetParameter(15,"fBkgMean",fParamStartValues[15], 1.e-08, 0., 1.e+04);
-               fitter->SetParameter(16,"fBkgSlope",fParamStartValues[16], 1.e-08, 0., 1.e+04);
-               fitter->SetParameter(17,"fBkgConst",fParamStartValues[17], 1.e-08, 0., 1.e+04);
+               fitter->SetParameter(14,"fBkgNorm",fParamStartValues[14], 1.e-08, -1.e+04, 1.e+04);
+               fitter->SetParameter(15,"fBkgMean",fParamStartValues[15], 1.e-08, -1.e+04, 1.e+04);
+               fitter->SetParameter(16,"fBkgSlope",fParamStartValues[16], 1.e-08, -1.e+04, 1.e+04);
+               fitter->SetParameter(17,"fBkgConst",fParamStartValues[17], 1.e-08, -1.e+04, 1.e+04);
                fitter->SetParameter(18,"fNormGaus1FF",fParamStartValues[18], 1.e-08, 0., 1.e+05);
                fitter->SetParameter(19,"fNormGaus2FF",fParamStartValues[19], 1.e-08, 0., 1.e+05); 
                fitter->SetParameter(20,"fMean1ResFF",fParamStartValues[20], 1.e-08, 0., 1.e+05);
@@ -241,15 +245,21 @@ Int_t AliDielectronBtoJPSItoEleCDFfitHandler::DoMinimization(Int_t step)
                 fitter->SetParameter(42,"fAlfaResSS",fParamStartValues[42], 1.e-08, 0., 1.e+05);
                 fitter->SetParameter(43,"fLambdaResSS",fParamStartValues[43], 1.e-08, 0., 1.e+05);
                 fitter->SetParameter(44,"fResNormExpSS",fParamStartValues[44], 1.e-08, 0., 1.e+05);
+                fitter->SetParameter(45,"fOneOvLamSym1",fParamStartValues[45], 1.e-10, 0.0000001, 5.e+01);
+                fitter->SetParameter(46,"fSym1",fParamStartValues[46], 1.e-08, 0., 1.e+06);
+                fitter->SetParameter(47,"fBkgInvMassPol4",fParamStartValues[47], 1.e-08, 0.0000001, 5.e+01);
+                fitter->SetParameter(48,"fBkgInvMassPol5",fParamStartValues[48], 1.e-08, 0., 5.e+06);
+                //(TMinuit*)fitter->GetMinuit()->SetErrorDef(fUp); 
                 }
 
-       for(UInt_t indexparam = 0; indexparam < 45; indexparam++){
+       for(UInt_t indexparam = 0; indexparam < 49; indexparam++){
                if(IsParamFixed(indexparam)) fitter->FixParameter((Int_t)indexparam); 
                else fitter->ReleaseParameter((Int_t)indexparam);
        }
-       Double_t arglist[2]={10000,0.1};
+       Double_t arglist[2]={10000,0.1}; Int_t iret = 0;
        if(step == 2) {Int_t  iret1 = fitter->ExecuteCommand("MINOS", arglist ,1); return iret1;}
-       Int_t iret=fitter->ExecuteCommand("MIGRAD", arglist ,2);
+       if(step == 0) { fitter->SetParameter(8,"fFsig",fParamStartValues[8], 1.e-10, 0., 1.); 
+                         iret=fitter->ExecuteCommand("MIGRAD", arglist ,2); }
        fitter->PrintResults(4,0);
 
        if(step == 3) {
@@ -288,7 +298,8 @@ Int_t AliDielectronBtoJPSItoEleCDFfitHandler::DoMinimization(Int_t step)
                 fContPlot2->Draw("l");
                 fContPlot1->Draw("l");
                
-                c2->Draw();    
+                c2->Draw();
+                c2->SaveAs("contourPlot.root");        
        }
 
        AliInfo("Minimization procedure finished\n");
@@ -308,7 +319,7 @@ void AliDielectronBtoJPSItoEleCDFfitHandler::CdfFCN(Int_t & /* npar */,
        TStopwatch t;
        t.Start();
 
-       f = fLikely->EvaluateLikelihood(fX,fM,fType,fNcand);
+       f = fLikely->EvaluateLikelihood(fX,fM,fPt, fType,fNcand);
 
        t.Stop();
        AliDebug(2,Form("Real time spent to calculate function == %f \n", t.RealTime()));
@@ -318,9 +329,9 @@ void AliDielectronBtoJPSItoEleCDFfitHandler::CdfFCN(Int_t & /* npar */,
        return;
 }
 //_______________________________________________________________________________________
-void AliDielectronBtoJPSItoEleCDFfitHandler::SetParamStartValues(Double_t inputparamvalues[45])
+void AliDielectronBtoJPSItoEleCDFfitHandler::SetParamStartValues(Double_t inputparamvalues[49])
 {
-       for(Int_t index=0; index < 45; index++) fParamStartValues[index] = inputparamvalues[index];
+       for(Int_t index=0; index < 49; index++) fParamStartValues[index] = inputparamvalues[index];
 }
 //_______________________________________________________________________________________
 void AliDielectronBtoJPSItoEleCDFfitHandler::SetResolutionConstants(Double_t* resolutionConst, Int_t type)
@@ -340,6 +351,18 @@ void AliDielectronBtoJPSItoEleCDFfitHandler::SetCrystalBallFunction(Bool_t okCB)
        //
        fLikely->SetCrystalBallFunction(okCB);
 }
+
+//_______________________________________________________________________________________
+void AliDielectronBtoJPSItoEleCDFfitHandler::SetExponentialFunction(Bool_t okExp)
+{
+        //
+        // Sets the CB as the parametrization for the signal invariant mass spectrum 
+        // (otherwise Landau is chosen)
+        //
+        fLikely->SetExponentialFunction(okExp);
+}
+
+
 //_______________________________________________________________________________________
 void AliDielectronBtoJPSItoEleCDFfitHandler::SetMassWndHigh(Double_t limit) 
 { 
index 4a0b710..787fcc8 100644 (file)
@@ -24,7 +24,7 @@ class AliDielectronBtoJPSItoEleCDFfitHandler : public TNamed {
                AliDielectronBtoJPSItoEleCDFfitHandler();\r
                AliDielectronBtoJPSItoEleCDFfitHandler& operator= (const  AliDielectronBtoJPSItoEleCDFfitHandler& c);\r
                AliDielectronBtoJPSItoEleCDFfitHandler(const AliDielectronBtoJPSItoEleCDFfitHandler& c);\r
-               AliDielectronBtoJPSItoEleCDFfitHandler(Double_t* decaytime, Double_t* invariantmass, Int_t *type, Int_t ncand);\r
+               AliDielectronBtoJPSItoEleCDFfitHandler(Double_t* decaytime, Double_t* invariantmass, Double_t *pt, Int_t *type, Int_t ncand);\r
                ~AliDielectronBtoJPSItoEleCDFfitHandler(); \r
                Double_t Up() const { return fUp; }\r
                void SetErrorDef(Double_t up) {fUp = up;}\r
@@ -38,10 +38,11 @@ class AliDielectronBtoJPSItoEleCDFfitHandler : public TNamed {
                 Double_t GetParameterError(Int_t numPar) const {return fitter->GetParError(numPar);}   \r
 \r
                 void FixParam(UInt_t param, Bool_t value) { fIsParamFixed.SetBitNumber(param,value); }\r
-               void FixAllParam(Bool_t value) { for(UInt_t par=0;par<45;par++) fIsParamFixed.SetBitNumber(par,value); }\r
+               void FixAllParam(Bool_t value) { for(UInt_t par=0;par<49;par++) fIsParamFixed.SetBitNumber(par,value); }\r
                Bool_t IsParamFixed(UInt_t param) { return fIsParamFixed.TestBitNumber(param); }\r
                void SetResolutionConstants(Double_t* resolutionConst, Int_t type);\r
                void SetCrystalBallFunction(Bool_t okCB);\r
+               void SetExponentialFunction(Bool_t okEpx);\r
                void SetMassWndHigh(Double_t limit);\r
                void SetMassWndLow(Double_t limit);\r
 \r
@@ -50,6 +51,7 @@ class AliDielectronBtoJPSItoEleCDFfitHandler : public TNamed {
 \r
                Double_t* Decaytime() const         { return fX; }\r
                Double_t* InvariantMass() const     { return fM; }\r
+               Double_t* TransverseMom() const     { return fPt; }\r
                 Int_t*    TypeCand() const          { return fType;}\r
                AliDielectronBtoJPSItoEleCDFfitFCN* LikelihoodPointer() const { return fLikely; }\r
                Int_t DoMinimization(Int_t step = 0);\r
@@ -59,10 +61,11 @@ class AliDielectronBtoJPSItoEleCDFfitHandler : public TNamed {
                //\r
                TBits fIsParamFixed;                               //array of bits: 0 = param free; 1 = param fixed;\r
                Bool_t fPrintStatus;                               //flag to enable the prit out of the algorithm at each step\r
-               Double_t fParamStartValues[45];                    //array of parameters input value\r
+               Double_t fParamStartValues[49];                    //array of parameters input value\r
                Double_t fUp;                                      //error definition \r
                Double_t* fX;                                      //pseudo-proper decay time X\r
                Double_t* fM;                                      //invariant mass M\r
+               Double_t* fPt;                                     //invariant mass M\r
                 Int_t* fType;                                      //candidate type\r
                AliDielectronBtoJPSItoEleCDFfitFCN* fLikely;       //Log likelihood function\r
                Int_t fNcand;                                      //number of candidates\r