Added D* (Raoul)
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 20 May 2011 22:07:27 +0000 (22:07 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 20 May 2011 22:07:27 +0000 (22:07 +0000)
PWG3/vertexingHF/AliAnalysisTaskSESignificance.cxx
PWG3/vertexingHF/AliAnalysisTaskSESignificance.h
PWG3/vertexingHF/AliHFMassFitter.cxx
PWG3/vertexingHF/macros/charmCutsOptimization.C

index 5775934..73a1144 100644 (file)
@@ -114,7 +114,12 @@ AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance(const char *name, T
 
   SetPDGCodes();
   SetDsChannel(kPhi);
-  SetMassLimits(0.15,fPDGmother); //check range
+  if (fDecChannel!=2) SetMassLimits(0.15,fPDGmother); //check range
+  else {
+    Float_t min = 0.13;
+    Float_t max = 0.19;
+    SetMassLimits(min, max);
+  }
   fNPtBins=fRDCuts->GetNPtBins();
 
   fNVars=fRDCuts->GetNVarsForOpt();
@@ -418,6 +423,7 @@ void AliAnalysisTaskSESignificance::UserCreateOutputObjects()
   fOutput->Add(fHistNEvents);
 
  
+  PostData(1,fOutput);    
   return;
 }
 
@@ -504,6 +510,18 @@ void AliAnalysisTaskSESignificance::UserExec(Option_t */*option*/)
   for (Int_t iProng = 0; iProng < nProng; iProng++) {
     fHistNEvents->Fill(6);
     d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
+
+    AliAODRecoCascadeHF* DStarToD0pi = NULL;
+    AliAODRecoDecayHF2Prong* D0Particle = NULL;
+    fPDGDStarToD0pi[0] = 421; fPDGDStarToD0pi[1] = 211;
+    fPDGD0ToKpi[0] = 321; fPDGD0ToKpi[1] = 211;
+
+    if (fDecChannel==2) {
+      DStarToD0pi = (AliAODRecoCascadeHF*)arrayProng->At(iProng);
+      if (!DStarToD0pi->GetSecondaryVtx()) continue;
+      D0Particle = (AliAODRecoDecayHF2Prong*)DStarToD0pi->Get2Prong();
+      if (!D0Particle) continue;
+      }
     
     Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(fPDGmother));
     Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel,aod);
@@ -564,7 +582,7 @@ void AliAnalysisTaskSESignificance::UserExec(Option_t */*option*/)
          FillD02p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
          break;
        case 2:
-         FillDstar(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
+         FillDstar(DStarToD0pi,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
          break;
        case 3:
          if(isSelected&1){
@@ -700,10 +718,56 @@ void AliAnalysisTaskSESignificance::FillD02p(AliAODRecoDecayHF* d,TClonesArray *
   }
 }
 
-void AliAnalysisTaskSESignificance::FillDstar(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Int_t /*matchtoMC*/){
+void AliAnalysisTaskSESignificance::FillDstar(AliAODRecoCascadeHF* dstarD0pi,TClonesArray *arrayMC,Int_t index,Int_t isSel){
     //D* channel
 
-  AliInfo("Dstar channel not implemented\n");
+    AliInfo("Dstar selected\n");
+    
+    Double_t mass = dstarD0pi->DeltaInvMass();
+
+    if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) fMassHist[index]->Fill(mass);
+    if(isSel>=2 && fPartOrAndAntiPart<=0) fMassHist[index]->Fill(mass);
+       
+    if(fReadMC) {
+       Int_t matchtoMC = -1; 
+       matchtoMC = dstarD0pi->MatchToMC(413,421,fPDGDStarToD0pi, fPDGD0ToKpi, arrayMC);
+
+       Int_t prongPdgDStarPlus=413,prongPdgDStarMinus=(-1)*413;
+       
+       if ((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) { 
+          //D*+
+         if(matchtoMC>=0) {
+             AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
+            Int_t pdgMC = dMC->GetPdgCode();
+       
+            if (pdgMC==prongPdgDStarPlus) fSigHist[index]->Fill(mass);
+            else {
+               dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
+               mass =  dstarD0pi->DeltaInvMass();
+               fRflHist[index]->Fill(mass);
+               dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
+             } 
+           } 
+         else fBkgHist[index]->Fill(mass);
+        }
+    
+       if (isSel>=2 && fPartOrAndAntiPart<=0) { 
+           //D*-
+          if (matchtoMC>=0) {
+             AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
+             Int_t pdgMC = dMC->GetPdgCode();
+       
+             if (pdgMC==prongPdgDStarMinus) fSigHist[index]->Fill(mass);
+             else {
+                dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
+                mass = dstarD0pi->DeltaInvMass();
+                fRflHist[index]->Fill(mass);
+                dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
+             }
+           } 
+           else fBkgHist[index]->Fill(mass);
+         }
+       }
 
 }
 
@@ -926,4 +990,3 @@ void AliAnalysisTaskSESignificance::Terminate(Option_t */*option*/)
   return;
 }
 //-------------------------------------------
-
index 2a94b71..a4df84f 100644 (file)
@@ -86,7 +86,7 @@ class AliAnalysisTaskSESignificance : public AliAnalysisTaskSE
   void FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel);
   void FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index, Int_t isSel);
   void FillDs(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel,Int_t optDecay);
-  void FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel);
+  void FillDstar(AliAODRecoCascadeHF* dstarD0pi,TClonesArray *arrayMC,Int_t index,Int_t isSel);
   void FillD04p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel);
   void FillLambdac(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index, Int_t isSel);
 
@@ -121,9 +121,10 @@ class AliAnalysisTaskSESignificance : public AliAnalysisTaskSE
   Int_t fNBins;  //number of bins in the mass histograms
   Int_t fPartOrAndAntiPart;  //fill histograms with particle only (+1), antiparticle only (-1), both (0)
   Int_t fDsChannel;          // Ds resonant channel selected
+  Int_t fPDGDStarToD0pi[2]; //PDG codes for the particles in the D* -> pi + D0 decay
+  Int_t fPDGD0ToKpi[2];    //PDG codes for the particles in the D0 -> K + pi decay
 
   ClassDef(AliAnalysisTaskSESignificance,4); // AliAnalysisTaskSE for the MC association of heavy-flavour decay candidates
 };
 
 #endif
-
index f95952d..d6b4ebe 100644 (file)
@@ -37,6 +37,7 @@
 #include <TMinuit.h>
 #include <TStyle.h>
 #include <TPaveText.h>
+#include <TDatabasePDG.h>
 
 #include "AliHFMassFitter.h"
 
@@ -246,6 +247,11 @@ void AliHFMassFitter::ComputeNFinalPars() {
     break;
   case 3:
     fNFinalPars=1;
+  case 4:
+    fNFinalPars=2;     
+    break;
+  case 5:
+    fNFinalPars=3;     
     break;
   default:
     cout<<"Error in computing fNFinalPars: check ftypeOfFit4Bkg"<<endl;
@@ -274,6 +280,12 @@ void AliHFMassFitter::ComputeParSize() {
   case 3:
     fParsSize = 1*3;
     break;
+  case 4:
+    fParsSize = 2*3; 
+    break;
+  case 5:
+    fParsSize = 3*3; 
+    break;
   default:
     cout<<"Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
     break;
@@ -608,6 +620,52 @@ Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){
     }
   }
 
+  //Power fit
+
+    // par[0] = tot integral
+    // par[1] = coef1
+    // par[2] = gaussian integral
+    // par[3] = gaussian mean
+    // par[4] = gaussian sigma
+
+  if (ftypeOfFit4Bkg==4) {
+    
+    if(ftypeOfFit4Sgn == 0) {
+      
+      Double_t parbkg[2] = {par[0]-par[2], par[1]}; 
+      bkg = FitFunction4Bkg(x,parbkg);
+    }
+    if(ftypeOfFit4Sgn == 1) {
+      
+      Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-par[2], par[1]}; 
+      bkg = FitFunction4Bkg(x,parbkg);
+    }
+    sgn = FitFunction4Sgn(x,&par[2]);
+  }
+
+
+  //Power and exponential fit
+
+    // par[0] = tot integral
+    // par[1] = coef1
+    // par[2] = coef2
+    // par[3] = gaussian integral
+    // par[4] = gaussian mean
+    // par[5] = gaussian sigma
+
+  if (ftypeOfFit4Bkg==5) {      
+   
+    if(ftypeOfFit4Sgn == 0) {
+       Double_t parbkg[3] = {par[0]-par[3],par[1],par[2]}; 
+       bkg = FitFunction4Bkg(x,parbkg); 
+    }
+    if(ftypeOfFit4Sgn == 1) {
+     Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-par[3], par[1], par[2]}; 
+     bkg = FitFunction4Bkg(x,parbkg);
+    }
+    sgn = FitFunction4Sgn(x,&par[3]);
+  }                            
+
   total = bkg + sgn;
   
   return  total;
@@ -683,12 +741,37 @@ Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){
   case 3:
     total=par[0+firstPar];
     break;
+  case 4:  
+    //power function 
+    //y=a(x-m_pi)^b -> integral = a/(b+1)*((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
+    //
+    //a = integral*(b+1)/((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
+    // * [0] = integralBkg;
+    // * [1] = b;
+    // a(power function) = [0]*([1]+1)/((max-m_pi)^([1]+1)-(min-m_pi)^([1]+1))*(x-m_pi)^[1]
+    {
+    Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
+
+    total = par[0+firstPar]*(par[1+firstPar]+1.)/(TMath::Power(fmaxMass-mpi,par[1+firstPar]+1.)-TMath::Power(fminMass-mpi,par[1+firstPar]+1.))*TMath::Power(x[0]-mpi,par[1+firstPar]);
+    }
+    break;
+  case 5:
+   //power function wit exponential
+    //y=a*Sqrt(x-m_pi)*exp(-b*(x-m_pi))  
+    { 
+    Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
+
+    total = par[1+firstPar]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[2+firstPar]*(x[0]-mpi));
+    } 
+    break;
 //   default:
 //     Types of Fit Functions for Background:
 //     * 0 = exponential;
 //     * 1 = linear;
 //     * 2 = polynomial 2nd order
 //     * 3 = no background"<<endl;
+//     * 4 = Power function 
+//     * 5 = Power function with exponential 
 
   }
   return total+gaus2;
@@ -937,6 +1020,14 @@ Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
       funcbkg->FixParameter(0,0.);
     }
     break;
+  case 4:
+    funcbkg->SetParNames("BkgInt","Coef2");  
+    funcbkg->SetParameters(sideBandsInt,0.5);
+    break;
+ case 5:
+    funcbkg->SetParNames("BkgInt","Coef1","Coef2");
+    funcbkg->SetParameters(sideBandsInt, 0.5, 50.);
+    break;
   default:
     cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
     return kFALSE;
@@ -963,6 +1054,9 @@ Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
     intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
     if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
+    if(ftypeOfFit4Bkg==5) conc1 = funcbkg->GetParameter(2); 
+
     //cout<<"First fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
   } 
   else cout<<"\t\t//"<<endl;
@@ -980,27 +1074,55 @@ Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
 
     funcbkg1->SetLineColor(2); //red
 
-    if(ftypeOfFit4Bkg==2){
-      cout<<"*** Polynomial Fit ***"<<endl;
-      funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
-      funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
-
-      //cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
-      //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
-    } else{
-      if(ftypeOfFit4Bkg==3) //no background: gaus sign+ gaus broadened
+ switch (ftypeOfFit4Bkg) {
+    case 0:
+       {
+        cout<<"*** Exponential Fit ***"<<endl;
+        funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
+       funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
+       }
+        break;
+     case 1: 
        {
-         cout<<"*** No background Fit ***"<<endl;
-         funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
-         funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.); 
-         funcbkg1->FixParameter(3,0.);
-       } else{ //expo or linear
-         if(ftypeOfFit4Bkg==0) cout<<"*** Exponential Fit ***"<<endl;
-         if(ftypeOfFit4Bkg==1) cout<<"*** Linear Fit ***"<<endl;
-         funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
-         funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
+       cout<<"*** Linear Fit ***"<<endl;
+        funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
+       funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
        }
+        break;    
+     case 2:
+        {
+        cout<<"*** Polynomial Fit ***"<<endl;
+        funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
+        funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
+        }
+        break;
+    case 3:
+       //no background: gaus sign+ gaus broadened
+       {
+       cout<<"*** No background Fit ***"<<endl;
+       funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
+       funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.); 
+       funcbkg1->FixParameter(3,0.);
+       } 
+        break;
+     case 4:
+       {       
+       cout<<"*** Power function Fit ***"<<endl;
+       funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef2");
+        funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
+               }
+        break;
+      case 5:
+       { 
+       cout<<"*** Power function conv. with exponential Fit ***"<<endl;
+        funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
+        funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
+       }
+        break;
     }
+      //cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
+      //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
+
     Int_t status=fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
     if (status != 0){
       cout<<"Minuit returned "<<status<<endl;
@@ -1017,6 +1139,8 @@ Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
     intbkg1=funcbkg1->GetParameter(3);
     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
     if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
+    if(ftypeOfFit4Bkg==5) conc1 = funcbkg1->GetParameter(5); 
+
 
   } else {
     bkgPar+=3;
@@ -1095,7 +1219,7 @@ Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
     if(fFixPar[2])funcmass->FixParameter(2,conc1);
     if(fFixPar[3])funcmass->FixParameter(3,sgnInt);
     if(fFixPar[4])funcmass->FixParameter(4,fMass);
-    if(fFixPar[5])funcmass->FixParameter(5,fSigmaSgn);
+    if(fFixPar[5])funcmass->FixParameter(5,fSigmaSgn); 
     //
     //funcmass->FixParameter(2,sgnInt);
   }
@@ -1241,13 +1365,21 @@ Bool_t AliHFMassFitter::RefitWithBkgOnly(Bool_t draw){
     funcbkg->SetParameters(integral,-10.,5);
     break;
   case 3:
-    cout<<"Warning! This choice does not have a lot of sense..."<<endl;
+    cout<<"Warning! This choice does not make a lot of sense..."<<endl;
     if(ftypeOfFit4Sgn==0){
       funcbkg->SetParNames("Const");
       funcbkg->SetParameter(0,0.);
       funcbkg->FixParameter(0,0.);
     }
     break;
+  case 4:     
+    funcbkg->SetParNames("BkgInt","Coef1");
+    funcbkg->SetParameters(integral,0.5);
+    break;
+  case 5:    
+    funcbkg->SetParNames("BkgInt","Coef1","Coef2");
+    funcbkg->SetParameters(integral,-10.,5.);
+    break;
   default:
     cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
     return kFALSE;
@@ -1269,6 +1401,7 @@ Bool_t AliHFMassFitter::RefitWithBkgOnly(Bool_t draw){
 }
 //_________________________________________________________________________
 Double_t AliHFMassFitter::GetChiSquare() const{
+  //Get Chi^2 method
   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
   if(!funcmass) {
     cout<<"funcmass not found"<<endl;
@@ -1279,6 +1412,7 @@ Double_t AliHFMassFitter::GetChiSquare() const{
 
 //_________________________________________________________________________
 Double_t AliHFMassFitter::GetReducedChiSquare() const{
+  //Get reduced Chi^2 method
   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
   if(!funcmass) {
     cout<<"funcmass not found"<<endl;
@@ -1499,6 +1633,12 @@ void AliHFMassFitter::WriteCanvas(TString userIDstring,TString path,Double_t nsi
   case 3:
     type="noB"; //3+1
     break;
+  case 4:  
+    type="Pow"; //3+3
+    break;
+  case 5:
+    type="PowExp"; //3+3
+    break;
   }
 
   TString filename=Form("%sMassFit.root",type.Data());
index 887a0b2..9745fd8 100644 (file)
@@ -1,6 +1,7 @@
 #include <fstream>
 #include <Riostream.h>
 #include <TSystem.h>
+#include <TMath.h>
 #include <TH1F.h>
 #include <TF1.h>
 #include <TFile.h>
 
 //global variables
 Double_t nsigma=3;
-Int_t decCh=0;
-Int_t fitbtype=0;
+Int_t decCh=2;
+Int_t fitbtype=5;
 Int_t rebin=2;
-Double_t sigma=0.012;
+Double_t sigma=0.0005;
 Int_t pdg;
 Double_t mass;
 
@@ -48,14 +49,22 @@ Bool_t MC(TH1F* hs,TH1F* hb, Double_t& sgn, Double_t& errsgn, Double_t& bkg, Dou
 //decCh:
 //- 0 = kDplustoKpipi
 //- 1 = kD0toKpi
-//- 2 = kDstartoKpipi
+//- 2 = kDStartoKpipi
 //- 3 = kDstoKKpi
 //- 4 = kD0toKpipipi
 //- 5 = kLambdactopKpi
 
 //Note: writefit=kTRUE writes the root files with the fit performed but it also draw all the canvas, so if your computer is not powerfull enough I suggest to run it in batch mode (root -b)
 
-Bool_t charmCutsOptimization(Bool_t isData=kTRUE,TString part="both"/*"A" anti-particle, "P" particle*/,TString centr="no",Bool_t writefit=kFALSE,Int_t minentries=50,Double_t *rangefit=0x0, Bool_t useBinCounting=kFALSE){
+// Method using this code:
+// root -b
+// .L macros/LoadLibraries
+// Loadlibraries()
+// .L macros/CharmCutsOptimization
+// charmCutsOptimization(/*options*/)
+
+Bool_t charmCutsOptimization(Bool_t isData=kTRUE,TString part="both"/*"A" anti-particle, "P" particle*/,TString centr="no",Bool_t writefit=kTRUE,Int_t minentries=50,Double_t *rangefit=0x0, Bool_t useBinCounting=kTRUE){
+
   outcheck.open("output.dat",ios::out);
   outdetail.open("discarddetails.dat",ios::out);
 
@@ -83,8 +92,8 @@ Bool_t charmCutsOptimization(Bool_t isData=kTRUE,TString part="both"/*"A" anti-p
     pdg=421;
     break;
   case 2:
-    listname+="Dstar";
-    mdvlistname+="Dstar";
+    listname+="Dstar0100";
+    mdvlistname+="Dstar0100";
     pdg=413;
     break;
   case 3:
@@ -107,6 +116,7 @@ Bool_t charmCutsOptimization(Bool_t isData=kTRUE,TString part="both"/*"A" anti-p
     return kFALSE;
   }
   mass=TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
+  if (decCh==2) mass=(TDatabasePDG::Instance()->GetParticle(pdg)->Mass() - TDatabasePDG::Instance()->GetParticle(421)->Mass()); 
 
   if(part!="both"){
     listname.Append(part);
@@ -293,8 +303,11 @@ Bool_t charmCutsOptimization(Bool_t isData=kTRUE,TString part="both"/*"A" anti-p
          cout<<name.Data()<<" not found"<<endl;
          continue;
        }
+
        MC(h,g,signal,errSignal,background,errBackground,signif,errSignif,sigmafit,status);
+
       }
+
       hstatus->Fill(status);
 
       if(status==1){
@@ -305,7 +318,6 @@ Bool_t charmCutsOptimization(Bool_t isData=kTRUE,TString part="both"/*"A" anti-p
        mdvB->SetElement(ih,background);
        mdvBerr->SetElement(ih,errBackground);
        hSigma[i]->Fill(sigmafit);
-      
       }else{
        mdvSgnf->SetElement(ih,0);
        mdvSgnferr->SetElement(ih,0);
@@ -322,18 +334,15 @@ Bool_t charmCutsOptimization(Bool_t isData=kTRUE,TString part="both"/*"A" anti-p
 
     }
   
-  
-
     fout->cd();
     mdvS->Write();
     mdvB->Write();
     mdvSgnf->Write();
-
     mdvSerr->Write();
     mdvBerr->Write();
     mdvSgnferr->Write();
     hSigma[i]->Write();
+
   }
   
 
@@ -341,17 +350,17 @@ Bool_t charmCutsOptimization(Bool_t isData=kTRUE,TString part="both"/*"A" anti-p
   cinfo->cd();
   cinfo->SetGrid();
   hstatus->Draw("htext0");
-
   fout->cd();
   hstatus->Write();
-
   fout->Close();
-
   outcheck<<"\nSummary:\n - Total number of histograms: "<<nhist<<"\n - "<<count<<" histograms with more than "<<minentries<<" entries; \n - Too few entries in histo "<<countNoHist<<" times;\n - Fit failed "<<countFitFail<<" times \n - no sense Signal/Background/Significance "<<countSgnfFail<<" times\n - only background "<<countBkgOnly<<" times"<<endl;
   outcheck.close();
+
+cout << "test1\n";
   return kTRUE;
 }
 
+
 //this function fit the hMass histograms
 //status = 0 -> fit fail
 //         1 -> fit ok and good results
@@ -364,7 +373,8 @@ Bool_t Data(TH1F* h,Double_t* rangefit,Bool_t writefit, Double_t& sgn, Double_t&
   Double_t min=h->GetBinLowEdge(7);
   Double_t max=h->GetBinLowEdge(nbin-5)+h->GetBinWidth(nbin-5);
 
-  min=h->GetBinLowEdge(1);
+  if(decCh != 2) min = h->GetBinLowEdge(1);
+  else min = TDatabasePDG::Instance()->GetParticle(211)->Mass();
   max=h->GetBinLowEdge(nbin+1);
 
   if(rangefit) {
@@ -490,14 +500,45 @@ Double_t ExpoBkgWoPeak(Double_t *x, Double_t *par){
 }
 
 //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+// par[0], par[1] expo params, par[2], par[3] exclusion range
+
+Double_t PowerBkgWoPeak(Double_t *x, Double_t *par){
+
+  if( reject && x[0]>par[2] && x[0]<par[3] ){
+    TF1::RejectPoint();
+    return 0;
+  }
+
+  Double_t xminusmpi = x[0]-TDatabasePDG::Instance()->GetParticle(211)->Mass();
+  return par[0]*TMath::Power(TMath::Abs(xminusmpi),par[1]) ;
+}
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+// par[0], par[1] expo params, par[2], par[3] exclusion range
+
+Double_t PowerExpoBkgWoPeak(Double_t *x, Double_t *par){
+
+  if( reject && x[0]>par[2] && x[0]<par[3] ){
+    TF1::RejectPoint();
+    return 0;
+  }
+  Double_t xminusmpi = x[0]-TDatabasePDG::Instance()->GetParticle(211)->Mass();
+  return par[0]*TMath::Sqrt(xminusmpi)*TMath::Exp(-1.*par[1]*xminusmpi);
+}
+
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 //this function fit the hMass histograms
 //status = 0 -> fit fail
 //         1 -> fit ok and good results
 //         2 -> negative signif
 Bool_t BinCounting(TH1F* h,Double_t* rangefit,Bool_t writefit, Double_t& sgn, Double_t& errsgn, Double_t& bkg, Double_t& errbkg, Double_t& sgnf, Double_t& errsgnf, Int_t& status){
 
+
+
   Int_t nbin=h->GetNbinsX();
-  Double_t min=h->GetBinLowEdge(1);
+  if(decCh != 2) Double_t min = h->GetBinLowEdge(1);
+  else Double_t min = TDatabasePDG::Instance()->GetParticle(211)->Mass();
   Double_t max=h->GetBinLowEdge(nbin+1);
 
   if(rangefit) {
@@ -505,9 +546,14 @@ Bool_t BinCounting(TH1F* h,Double_t* rangefit,Bool_t writefit, Double_t& sgn, Do
     max=rangefit[1];
   }
 
-  // Bkg fit : exponential = A*exp(B*x) 
   reject = true;
-  TF1 *fBkgFit = new TF1("fBkgFit",ExpoBkgWoPeak,min,max,4);
+//Bkg fit : exponential = A*exp(B*x) 
+  if(decCh != 2) TF1 *fBkgFit = new TF1("fBkgFit",ExpoBkgWoPeak,min,max,4);          
+  else { 
+   if (fitbtype == 5) TF1 *fBkgFit = new TF1("fBkgFit",PowerExpoBkgWoPeak,min,max,4); //Bkg fit : PowerExpo = A*(x-mpi)^(1/2)*exp(-B*(x-mpi))
+   else TF1 *fBkgFit = new TF1("fBkgFit",PowerBkgWoPeak,min,max,4); //Bkg fit : Power = A*(x-mpi)^B
+  }
+
   fBkgFit->FixParameter(2,mass-nsigma*sigma);
   fBkgFit->FixParameter(3,mass+nsigma*sigma);
   TFitResultPtr r = h->Fit(fBkgFit,"RS+");
@@ -521,12 +567,17 @@ Bool_t BinCounting(TH1F* h,Double_t* rangefit,Bool_t writefit, Double_t& sgn, Do
   } 
 
   reject = false;
-  TF1 *fBkgFct = new TF1("fBkgFct",ExpoBkgWoPeak,min,max,4);
+  if(decCh !=2) TF1 *fBkgFct = new TF1("fBkgFct",ExpoBkgWoPeak,min,max,4);  
+  else {
+   if (fitbtype == 5) TF1 *fBkgFct = new TF1("fBkgFct",PowerExpoBkgWoPeak,min,max,4);  
+   else TF1 *fBkgFct = new TF1("fBkgFct",PowerBkgWoPeak,min,max,4);
+  }
   fBkgFct->SetLineStyle(2);
   for(Int_t i=0; i<4; i++) fBkgFct->SetParameter(i,fBkgFit->GetParameter(i));
   h->GetListOfFunctions()->Add(fBkgFct);
   TH1F * hBkgFct = (TH1F*)fBkgFct->GetHistogram();
 
+
   status = 1;    
   Double_t binStartCount = h->FindBin(mass-nsigma*sigma);
   Double_t binEndCount = h->FindBin(mass+nsigma*sigma);
@@ -535,9 +586,11 @@ Bool_t BinCounting(TH1F* h,Double_t* rangefit,Bool_t writefit, Double_t& sgn, Do
     counts += h->GetBinContent( ibin );
     errcounts += counts ;
     Double_t center =  h->GetBinCenter(ibin);
-    bkgcounts += hBkgFct->GetBinContent( hBkgFct->FindBin(center) );
+    bkgcounts += hBkgFct->GetBinContent( hBkgFct->FindBin(center) ); 
+    
     errbkgcounts += bkgcounts ;
   }
+  
   bkg = bkgcounts;
   errbkg = TMath::Sqrt( errbkgcounts );
   sgn = counts - bkg ;
@@ -545,14 +598,20 @@ Bool_t BinCounting(TH1F* h,Double_t* rangefit,Bool_t writefit, Double_t& sgn, Do
   errsgn = TMath::Sqrt( counts + errbkg*errbkg );
   sgnf = sgn / TMath::Sqrt( sgn + bkg );
   errsgnf = TMath::Sqrt( sgnf*sgnf/(sgn+bkg)/(sgn+bkg)*(1/4.*errsgn*errsgn+errbkg*errbkg) + sgnf*sgnf/sgn/sgn*errsgn*errsgn );
-  //    cout << " Signal "<<sgn<<" +- "<<errsgn<<", bkg "<<bkg<<" +- "<<errbkg<<", significance "<<sgnf<<" +- "<<errsgnf<<endl;
+      cout << " Signal "<<sgn<<" +- "<<errsgn<<", bkg "<<bkg<<" +- "<<errbkg<<", significance "<<sgnf<<" +- "<<errsgnf<<endl;
   
   if(writefit) {
+
     TString filename = Form("%sMassFit.root",h->GetName());
+
     TFile* outputcv = new TFile(filename.Data(),"recreate");      
+
     TCanvas* c = new TCanvas();
+
     c->SetName(Form("%s",h->GetName()));
+
     h->Draw();
+
     TPaveText *pavetext=new TPaveText(0.4,0.7,0.85,0.9,"NDC");     
     pavetext->SetBorderSize(0);
     pavetext->SetFillStyle(0);
@@ -560,18 +619,27 @@ Bool_t BinCounting(TH1F* h,Double_t* rangefit,Bool_t writefit, Double_t& sgn, Do
     pavetext->AddText(Form("Bkg = %4.2e #pm %4.2e",bkg,errbkg));
     pavetext->AddText(Form("Signif = %3.2f #pm %3.2f",sgnf,errsgnf));
     c->cd();
+
     pavetext->DrawClone();
+
     outputcv->cd();
+
     c->Write();
+
+
     outputcv->Close();
+
     delete outputcv;
+
     delete c;
   }
 
-  
+
   delete fBkgFit;
+
   delete fBkgFct;
 
+
   return kTRUE; 
 }
 
@@ -1030,7 +1098,7 @@ Float_t* GetCutValuesFromNHist(AliMultiDimVector* vct,Int_t ptbin,Int_t nhist){
 // valsgiven = array of dimention = to the number of variables optimized. For each variable put kTRUE if the value is given (in values), kFALSE otherwise
 // 
 
-void DrawPossibilities(Int_t ptbin,Bool_t* valsgiven,Float_t* values,TString path="./",Int_t decCh=1){
+void DrawPossibilities(Int_t ptbin,Bool_t* valsgiven,Float_t* values,TString path="./",Int_t decCh=2){
   gStyle->SetFrameBorderMode(0);
   gStyle->SetCanvasColor(0);
   gStyle->SetFrameFillColor(0);
@@ -1063,8 +1131,8 @@ void DrawPossibilities(Int_t ptbin,Bool_t* valsgiven,Float_t* values,TString pat
     
     break;
   case 2:
-    listname+="Dstar";
-    mdvlistname+="Dstar";
+    listname+="Dstar0100";
+    mdvlistname+="Dstar0100";
     
     break;
   case 3:
@@ -1121,7 +1189,7 @@ void DrawPossibilities(Int_t ptbin,Bool_t* valsgiven,Float_t* values,TString pat
   }
 }
 
-void Merge2Bins(Int_t b1, Int_t b2,TString pathin="./",Int_t decCh=1,TString part="both"/*"A" anti-particle, "P" particle*/){
+void Merge2Bins(Int_t b1, Int_t b2,TString pathin="./",Int_t decCh=2,TString part="both"/*"A" anti-particle, "P" particle*/){
 
   if(b2!=b1+1){
     printf("The bins to be merget must be consecutive. Check! [b1 = %d, b2= %d]\n",b1,b2);
@@ -1146,8 +1214,8 @@ void Merge2Bins(Int_t b1, Int_t b2,TString pathin="./",Int_t decCh=1,TString par
     mdvlistname+="D0";
     break;
   case 2:
-    listname+="Dstar";
-    mdvlistname+="Dstar";
+    listname+="Dstar0100";
+    mdvlistname+="Dstar0100";
     break;
   case 3:
     listname+="Ds";
@@ -1392,3 +1460,4 @@ void SubtractBkg(Int_t nhisto){
   c->SaveAs(Form("h%d%sSubtr.png",nhisto,fitType.Data()));
   cvnewfits->SaveAs(Form("h%d%sFitNew.png",nhisto,fitType.Data()));
 }
+