Update of B->Jpsi->ee analysis classes
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 12 Oct 2009 12:39:01 +0000 (12:39 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 12 Oct 2009 12:39:01 +0000 (12:39 +0000)
PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx
PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h
PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx
PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h
PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx
PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h
PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx
PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h
PWG3/vertexingHF/macros/FitCDFLocal.C [new file with mode: 0644]

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