Signal extraction and error calc. updates (Ionut), updated configs (Jens)
authorandronic <andronic@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Sep 2010 17:41:20 +0000 (17:41 +0000)
committerandronic <andronic@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Sep 2010 17:41:20 +0000 (17:41 +0000)
PWG3/dielectron/AliDielectronSignalBase.cxx
PWG3/dielectron/AliDielectronSignalBase.h
PWG3/dielectron/AliDielectronSignalExt.cxx
PWG3/dielectron/AliDielectronSignalExt.h
PWG3/dielectron/AliDielectronSignalFunc.cxx
PWG3/dielectron/AliDielectronSignalFunc.h
PWG3/dielectron/AliDielectronSpectrum.cxx
PWG3/dielectron/macros/ConfigJpsi2eeData.C
PWG3/dielectron/macros/ConfigJpsi2eeFilter.C

index da5386b..f4c5157 100644 (file)
@@ -30,29 +30,49 @@ resulting from single and mixed events, as defined in AliDielectron.cxx
 
 #include "TPaveText.h"
 #include "AliDielectronSignalBase.h"
+#include <TH1F.h>
 
 ClassImp(AliDielectronSignalBase)
 
 AliDielectronSignalBase::AliDielectronSignalBase() :
   TNamed(),
+  fHistSignal(0),
+  fHistBackground(0),
+  fHistDataPM(0),
+  fHistDataPP(0),
+  fHistDataMM(0),
   fValues(6),
   fErrors(6),
-  fIntMin(3.01),
-  fIntMax(3.19)
+  fIntMin(0),
+  fIntMax(0),
+  fFitMin(0),
+  fFitMax(0),
+  fRebin(1),
+  fMethod(kLikeSign),
+  fProcessed(kFALSE)
 {
   //
   // Default Constructor
   //
-  
 }
 
 //______________________________________________
 AliDielectronSignalBase::AliDielectronSignalBase(const char* name, const char* title) :
   TNamed(name, title),
+  fHistSignal(0),
+  fHistBackground(0),
+  fHistDataPM(0),
+  fHistDataPP(0),
+  fHistDataMM(0),
   fValues(6),
   fErrors(6),
-  fIntMin(3.01),
-  fIntMax(3.19)
+  fIntMin(0),
+  fIntMax(0),
+  fFitMin(0),
+  fFitMax(0),
+  fRebin(1),
+  fMethod(kLikeSign),
+  fProcessed(kFALSE)
 {
   //
   // Named Constructor
@@ -65,7 +85,8 @@ AliDielectronSignalBase::~AliDielectronSignalBase()
   //
   // Default Destructor
   //
-  
+  if(fHistSignal) delete fHistSignal;
+  if(fHistBackground) delete fHistBackground;
 }
 
 //______________________________________________
@@ -85,13 +106,14 @@ TPaveText* AliDielectronSignalBase::DrawStats(Double_t x1/*=0.*/, Double_t y1/*=
   t->SetFillColor(kWhite);
   t->SetBorderSize(1);
   t->SetTextAlign(12);
-  t->AddText(Form("Signal : %.5g #pm %.5g",GetSignal(),GetSignalError()));
-  t->AddText(Form("Backgnd: %.5g #pm %.5g",GetBackground(),GetBackgroundError()));
-  t->AddText(Form("Signif.: %.5g #pm %.5g",GetSignificance(),GetSignificanceError()));
-  t->AddText(Form("SoB    : %.5g #pm %.5g",GetSignalOverBackground(),GetSignalOverBackgroundError()));
-  if (GetMass()>0){
-    t->AddText(Form("Mass: %.5g #pm %.5g", GetMass(), GetMassError()));
-    t->AddText(Form("Mass res.: %.5g #pm %.5g", GetMassWidth(), GetMassWidthError()));
+  t->AddText(Form("Range  : %.2f - %.2f GeV/c^{2}", fIntMin, fIntMax));
+  t->AddText(Form("Signal : %.1f #pm %.1f", fValues(0), fErrors(0)));
+  t->AddText(Form("Backgnd: %.1f #pm %.1f", fValues(1), fErrors(1)));
+  t->AddText(Form("Signif.: %.2f #pm %.2f", fValues(2), fErrors(2)));
+  t->AddText(Form("S/B    : %.2f #pm %.2f", fValues(3), fErrors(3)));
+  if(fValues(4)>0) {
+    t->AddText(Form("Mass: %.2f #pm %.2f GeV/c^{2}", fValues(4), fErrors(4)));
+    t->AddText(Form("Mass res.: %.1f #pm %.1f MeV/c^{2}", 1000*fValues(5), 1000*fErrors(5)));
   }
   t->Draw();
 
@@ -104,16 +126,12 @@ void AliDielectronSignalBase::Print(Option_t */*option*/) const
   //
   // Print the statistics
   //
-  printf("Signal : %.5g #pm %.5g\n",GetSignal(),GetSignalError());
-  printf("Backgnd: %.5g #pm %.5g\n",GetBackground(),GetBackgroundError());
-  printf("Signif.: %.5g #pm %.5g\n",GetSignificance(),GetSignificanceError());
-  printf("SoB    : %.5g #pm %.5g\n",GetSignalOverBackground(),GetSignalOverBackgroundError());
-  if (GetMass()>0){
-    printf("Mass: %.5g #pm %.5g\n", GetMass(), GetMassError());
-    printf("Mass res.: %.5g #pm %.5g\n", GetMassWidth(), GetMassWidthError());
+  printf("Signal : %.5g #pm %.5g\n",fValues(0), fErrors(0));
+  printf("Backgnd: %.5g #pm %.5g\n",fValues(1), fErrors(1));
+  printf("Signif.: %.5g #pm %.5g\n",fValues(2), fErrors(2));
+  printf("SoB    : %.5g #pm %.5g\n",fValues(3), fErrors(3));
+  if(fValues(4)>0){
+    printf("Mass: %.5g #pm %.5g\n", fValues(4), fErrors(4));
+    printf("Mass res.: %.5g #pm %.5g\n", fValues(5), fErrors(5));
   }
 }
-
-
-
-
index 139a5d1..5be8247 100644 (file)
 //#                                                           #
 //#############################################################
 
+
 #include <TNamed.h>
 #include <TVectorT.h>
 #include <TMath.h>
 
 class TObjArray;
 class TPaveText;
-class TH1;
+class TH1F;
 
 class AliDielectronSignalBase : public TNamed {
 public:
+  enum EBackgroundMethod {
+    kFitted = 0,
+    kLikeSign,
+    kEventMixing
+  };
+
   AliDielectronSignalBase();
   AliDielectronSignalBase(const char*name, const char* title);
-
+  
   virtual ~AliDielectronSignalBase();
-
   
   void SetIntegralRange(Double_t min, Double_t max) {fIntMin=min;fIntMax=max;}
-  
-  void SetSignal(Double_t val, Double_t valErr)               {fValues(0)=val; fErrors(0)=valErr;}
-  void SetBackground(Double_t val, Double_t valErr)           {fValues(1)=val; fErrors(1)=valErr;}
-  void SetSignificance(Double_t val, Double_t valErr)         {fValues(2)=val; fErrors(2)=valErr;}
-  void SetSignalOverBackground(Double_t val, Double_t valErr) {fValues(3)=val; fErrors(3)=valErr;}
-  void SetMass(Double_t val, Double_t valErr)                 {fValues(4)=val; fErrors(4)=valErr;}
-  void SetMassWidth(Double_t val, Double_t valErr)            {fValues(5)=val; fErrors(5)=valErr;}
-  
+  void SetFitRange(Double_t min, Double_t max) {fFitMin=min; fFitMax=max;}
+  void SetRebin(Int_t factor) {fRebin = factor;}
+  void SetMethod(EBackgroundMethod method) {fMethod = method;}
+
   const TVectorD& GetValues() const {return fValues;}
   const TVectorD& GetErrors() const {return fErrors;}
 
-  Double_t GetIntegralMin() const { return fIntMin; }
-  Double_t GetIntegralMax() const { return fIntMax; }
-  
-  Double_t GetSignal()               const {return fValues(0);}
-  Double_t GetBackground()           const {return fValues(1);}
-  Double_t GetSignificance()         const {return fValues(2);}
-  Double_t GetSignalOverBackground() const {return fValues(3);}
-  Double_t GetMass()                 const {return fValues(4);}
-  Double_t GetMassWidth()            const {return fValues(5);}
-  
-  Double_t GetSignalError()               const {return fErrors(0);}
-  Double_t GetBackgroundError()           const {return fErrors(1);}
-  Double_t GetSignificanceError()         const {return fErrors(2);}
-  Double_t GetSignalOverBackgroundError() const {return fErrors(3);}
-  Double_t GetMassError()                 const {return fErrors(4);}
-  Double_t GetMassWidthError()            const {return fValues(5);}
-  
-  void GetSignal(Double_t &val, Double_t &valErr)               {val=fValues(0); valErr=fErrors(0);}
-  void GetBackground(Double_t &val, Double_t &valErr)           {val=fValues(1); valErr=fErrors(1);}
-  void GetSignificance(Double_t &val, Double_t &valErr)         {val=fValues(2); valErr=fErrors(2);}
-  void GetSignalOverBackground(Double_t &val, Double_t &valErr) {val=fValues(3); valErr=fErrors(3);}
+  Double_t GetIntegralMin()          const { return fIntMin; }
+  Double_t GetIntegralMax()          const { return fIntMax; }
+  Double_t GetSignal()               const { return fValues(0);}
+  Double_t GetSignalError()          const { return fErrors(0);}
+  Double_t GetBackground()           const { return fValues(1);}
+  Double_t GetBackgroundError()      const { return fErrors(1);}
+  Double_t GetSignificance()         const { return fValues(2);}
+  Double_t GetSignificanceError()    const { return fErrors(2);}
+  Double_t GetSB()                   const { return fValues(3);}
+  Double_t GetSBError()              const { return fErrors(3);}
+  Double_t GetMass()                 const { return fValues(4);}
+  Double_t GetMassError()            const { return fErrors(4);}
+  Double_t GetMassWidth()            const { return fValues(5);}
+  Double_t GetMassWidthError()       const { return fErrors(5);}
+
+  TH1F* GetSignalHistogram()        {return fHistSignal;}
+  TH1F* GetBackgroundHistogram()    {return fHistBackground;}
+
+  virtual void Print(Option_t *option="") const;
 
   /**
   This function needs to be implemented by the signal extraction classes.
+  Here all the work should be done.
   
   The signal extraction is done on the mass spectra.
   The TObjArray should contain the Inv. Mass spectra of the 10 possible combinations
      for single and mixed events defined in AliDielectron.cxx
-  In the referece TVectorDs the values and value errors of the result should be stored:
-  Parameter - 0: Signal entries
-              1: Background entries
-              2: Significance ( S/sqr(S+B) )
-              3: Signal/Background
-              4: Mass
-              5: Mass width
-
-  It is enough to calculate the signal and background and then call
-            SetSignificanceAndSOB(TVectorD &values, TVectorD &errors)
-            Please also read the description there!!!
   */
   virtual void Process(TObjArray * const /*arrhist*/) = 0;
+    
+protected: 
 
-  virtual void Print(Option_t *option="") const;
-  
-protected:
-  
-  void SetSignificanceAndSOB();
+  TH1F *fHistSignal;                  // histogram of pure signal
+  TH1F *fHistBackground;              // histogram of background (fitted=0, like-sign=1, event mixing=2)
+  TH1F *fHistDataPM;                  // histogram of selected +- pair candidates
+  TH1F *fHistDataPP;                  // histogram of selected ++ pair candidates
+  TH1F *fHistDataMM;                  // histogram of selected -- pair candidates
+
+  TVectorD fValues;                   // values
+  TVectorD fErrors;                   // value errors
+
+  Double_t fIntMin;                   // signal extraction range min
+  Double_t fIntMax;                   // signal extraction range max
+  Double_t fFitMin;                   // fit range lowest inv. mass
+  Double_t fFitMax;                   // fit range highest inv. mass
+
+  Int_t fRebin;                       // histogram rebin factor
+  EBackgroundMethod fMethod;          // method for background substraction
+
+  Bool_t fProcessed;                  // flag
+
+  void SetSignificanceAndSOB();       // calculate the significance and S/B
   TPaveText* DrawStats(Double_t x1=0., Double_t y1=0., Double_t x2=0., Double_t y2=0.);
-  void Reset() {fValues.Zero(); fErrors.Zero();}
-  
-private:
-  TVectorD fValues;            // values
-  TVectorD fErrors;            // value errors
-  
-  Double_t fIntMin;            // signal extraction range min
-  Double_t fIntMax;            // signal extraction range max
-  
+
   AliDielectronSignalBase(const AliDielectronSignalBase &c);
   AliDielectronSignalBase &operator=(const AliDielectronSignalBase &c);
 
-  ClassDef(AliDielectronSignalBase,2)         // Dielectron SignalBase
+  ClassDef(AliDielectronSignalBase,3) // Dielectron SignalBase
 };
 
-//
-// Inline functions
-//
-
 inline void AliDielectronSignalBase::SetSignificanceAndSOB()
 {
   //
-  // calculate significance and signal over background from values
-  // it is expected that:
-  // - 'values' and 'errors' have the size 4
-  // - Parameters 0 are given: signal and signal error
-  // - Parameters 1 are given: background and background error
-  // - Optionally parameter 2 can be given: signal+background and its error
+  // Calculate S/B and significance
   //
-  // The calculated significance and S/B and the corresponding errors
-  //   are written in parameters 2 and 3, the first two parameters (signal and background)
-  //   stay untouched
-
-  Double_t signal=fValues(0),      signal_err=fErrors(0);
-  Double_t background=fValues(1),  background_err=fErrors(1);
-  Double_t sigPlusBack=fValues(2), sigPlusBack_err=fErrors(2);
-  Double_t significance=0,        significance_err=0;
-  Double_t sob=0,                 sob_err=0;
-  
-  if (sigPlusBack<1e-20){
-    //calculate signal plus background
-    sigPlusBack=signal+background;
-    sigPlusBack_err=TMath::Sqrt( signal_err*signal_err + background_err*background_err );
-  }
-
-  if (sigPlusBack>1e-30){
-
-    significance=signal/TMath::Sqrt(sigPlusBack);
-    
-    significance_err=TMath::Sqrt((signal_err/signal)*(signal_err/signal)+
-                                 (0.5*sigPlusBack_err/sigPlusBack)*(0.5*sigPlusBack_err/sigPlusBack)
-                                )*significance;
-    
-  }
-  
-  if (background>1e-30){
-    sob=signal/background;
-    sob_err=TMath::Sqrt((signal_err/signal)*(signal_err/signal)+
-                        (background_err/background)*(background_err/background))*sob;
-  }
-
-  fValues(2)=significance;
-  fErrors(2)=significance_err;
-  
-  fValues(3)=sob;
-  fErrors(3)=sob_err;
+  // Signal/Background
+  fValues(3) = (fValues(1)>0 ? fValues(0)/fValues(1) : 0);
+  Float_t epsSig = (fValues(0)>0 ? fErrors(0)/fValues(0) : 0);
+  Float_t epsBknd = (fValues(1)>0 ? fErrors(1)/fValues(1) : 0);
+  fErrors(3) = fValues(3)*TMath::Sqrt(epsSig*epsSig + epsBknd*epsBknd);
+  // Significance
+  fValues(2) = ((fValues(0)+fValues(1))>0 ? fValues(0)/TMath::Sqrt(fValues(0)+fValues(1)) : 0);
+  Float_t s = fValues(0); Float_t b = fValues(1);
+  fErrors(2) = ((s+b)>0 ? TMath::Sqrt((s*(s+2*b)*(s+2*b)+b*s*s)/(4*TMath::Power(s+b,3))) : 0);
 }
 
 #endif
index 0d12f9b..75b6aef 100644 (file)
@@ -27,12 +27,15 @@ can be used.
 */
 //                                                                       //
 ///////////////////////////////////////////////////////////////////////////
-
 #include <TF1.h>
 #include <TH1.h>
+#include <TH2F.h>
+#include <TLatex.h>
+#include <TLegend.h>
 #include <TCanvas.h>
 #include <TMath.h>
 #include <TString.h>
+#include <TLine.h>
 
 #include <AliLog.h>
 
@@ -41,41 +44,16 @@ can be used.
 ClassImp(AliDielectronSignalExt)
 
 AliDielectronSignalExt::AliDielectronSignalExt() :
-  AliDielectronSignalBase(),
-  fSignPM(0x0),
-  fSignPP(0x0),
-  fSignMM(0x0),
-  fBackground(0x0),
-  fSignal(0x0),
-  fMethod(1),
-  fRebin(1),
-  fBins(0),
-  fDrawMin(0.),
-  fDrawMax(0.),
-  fFitMin(2.5),
-  fFitMax(4)
+  AliDielectronSignalBase()
 {
   //
   // Default Constructor
   //
-  
 }
 
 //______________________________________________
 AliDielectronSignalExt::AliDielectronSignalExt(const char* name, const char* title) :
-  AliDielectronSignalBase(name, title),
-  fSignPM(0x0),
-  fSignPP(0x0),
-  fSignMM(0x0),
-  fBackground(0x0),
-  fSignal(0x0),
-  fMethod(1),
-  fRebin(1),
-  fBins(0),
-  fDrawMin(0.),
-  fDrawMax(0.),
-  fFitMin(2.5),
-  fFitMax(4)
+  AliDielectronSignalBase(name, title)
 {
   //
   // Named Constructor
@@ -88,22 +66,6 @@ AliDielectronSignalExt::~AliDielectronSignalExt()
   //
   // Default Destructor
   //
-  if (fSignPM)     delete fSignPM;
-  if (fSignPP)     delete fSignPP;
-  if (fSignMM)     delete fSignMM;
-  if (fBackground) delete fBackground;
-  if (fSignal)     delete fSignal;
-}
-
-//______________________________________________
-void AliDielectronSignalExt::SetHistograms(TH1F* const unlike, TH1F* const backg, TH1F* const signal)
-{
-  //
-  // set histograms 
-  //
-  fSignPM = (TH1F*)unlike->Clone("fSignPM");
-  fBackground = backg;
-  fSignal = signal;
 }
 
 //______________________________________________
@@ -113,11 +75,11 @@ void AliDielectronSignalExt::Process(TObjArray* const arrhist)
   // signal subtraction. support like-sign subtraction and event mixing method
   //
   switch ( fMethod ){
-    case 1 :
+    case kLikeSign :
       ProcessLS(arrhist);    // process like-sign subtraction method
       break;
 
-    case 2 : 
+    case kEventMixing : 
       ProcessEM(arrhist);    // process event mixing method
       break;
 
@@ -132,53 +94,55 @@ void AliDielectronSignalExt::ProcessLS(TObjArray* const arrhist)
   //
   // signal subtraction 
   //
-  fSignPP = (TH1F*)arrhist->At(0);  // like sign   : plus-plus
-  fSignPM = (TH1F*)arrhist->At(1);  // unlike-sign : plus-minus
-  fSignMM = (TH1F*)arrhist->At(2);  // like sign   : minus-minus
-  fSignPP->Sumw2();
-  fSignPM->Sumw2();
-  fSignMM->Sumw2();
+  fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP");  // ++    SE
+  fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
+  fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM");  // --    SE   
+  fHistDataPP->Sumw2();
+  fHistDataPM->Sumw2();
+  fHistDataMM->Sumw2();
   
-  if ( fRebin>1 ){ Rebin(fRebin); }       // rebinning of histogram
-
-  fBins = fSignPM->GetNbinsX();            // number of bins
-  Double_t minX  = fSignPM->GetBinLowEdge(1);       // minimum X value in axis
-  Double_t maxX  = fSignPM->GetBinLowEdge(fBins+1); // maximum X value in axis
-  
-  AliDebug(10,Form("histogram #bin = %d , min = %f , max = %f\n",fBins,minX,maxX));
-  TH1F* hBackground = new TH1F("hBackground","Like-sign background",fBins,minX,maxX); 
-  TH1F* hSubtracted = new TH1F("hSubtracted","Background subtracted",fBins,minX,maxX);
+  // rebin the histograms
+  if (fRebin>1) { 
+    fHistDataPP->Rebin(fRebin);
+    fHistDataPM->Rebin(fRebin);
+    fHistDataMM->Rebin(fRebin);
+  }       
+
+  fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal", 
+                        fHistDataPM->GetXaxis()->GetNbins(),
+                        fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
+  fHistBackground = new TH1F("HistBackground", "Like-sign contribution", 
+                            fHistDataPM->GetXaxis()->GetNbins(),
+                            fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 
   // fill out background and subtracted histogram
-  for ( Int_t ibin=1; ibin<fBins+1; ibin++ ){
-    Double_t mass  = ibin*(maxX-minX)/(Double_t)fBins;
-    Double_t backgr = 2.*sqrt( fSignPP->GetBinContent(ibin) * fSignMM->GetBinContent(ibin) );
-    Double_t signal = fSignPM->GetBinContent(ibin) - backgr;
-
-    hBackground->Fill(mass, backgr);
-    hSubtracted->Fill(mass, signal);
+  for(Int_t ibin=1; ibin<=fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
+    Float_t pm = fHistDataPM->GetBinContent(ibin);
+    Float_t pp = fHistDataPP->GetBinContent(ibin);
+    Float_t mm = fHistDataMM->GetBinContent(ibin);
+    Float_t epm = fHistDataPM->GetBinError(ibin);
+
+    Float_t background = 2*TMath::Sqrt(pp*mm);
+    Float_t ebackground = TMath::Sqrt(mm+pp);
+    Float_t signal = pm - background;
+    Float_t error = TMath::Sqrt(epm*epm+mm+pp);
+
+    fHistSignal->SetBinContent(ibin, signal);
+    fHistSignal->SetBinError(ibin, error);
+    fHistBackground->SetBinContent(ibin, background);
+    fHistBackground->SetBinError(ibin, ebackground);
   }
-  SetHistograms(fSignPM, hBackground, hSubtracted);
-
-
-  Double_t signal=0, signal_err=0;
-  Double_t background=0, background_err=0;
-  
-  signal     = fSignal->IntegralAndError(fSignal->FindBin(GetIntegralMin()),
-                                         fSignal->FindBin(GetIntegralMax()), signal_err);
-  background = fBackground->IntegralAndError(fBackground->FindBin(GetIntegralMin()),
-                                             fBackground->FindBin(GetIntegralMax()), background_err);
-
-  //reset result arrays
-  Reset();
-  //set values
-  SetSignal(signal,signal_err);
-  SetBackground(background,background_err);
+  // signal
+  fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
+                                            fHistSignal->FindBin(fIntMax), fErrors(0));
+  // background
+  fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
+                                                 fHistBackground->FindBin(fIntMax), 
+                                                 fErrors(1));
+  // S/B and significance
   SetSignificanceAndSOB();
 
-  // cleanup
-  //delete hBackground;
-  //delete hSubtracted;
+  fProcessed = kTRUE;
 }
 
 //______________________________________________
@@ -192,17 +156,6 @@ void AliDielectronSignalExt::ProcessEM(TObjArray* const arrhist)
 }
 
 //______________________________________________
-void AliDielectronSignalExt::Rebin(Int_t rebin)
-{
-  // 
-  // rebinning of histograms
-  //
-  fSignPM->Rebin(rebin);
-  fSignPP->Rebin(rebin);
-  fSignMM->Rebin(rebin);
-}
-
-//______________________________________________
 void AliDielectronSignalExt::Draw(const Option_t* option)
 {
   //
@@ -211,38 +164,124 @@ void AliDielectronSignalExt::Draw(const Option_t* option)
   TString drawOpt(option); 
   drawOpt.ToLower();   
 
+  Float_t minY = 0.001;
+  Float_t maxY = 1.2*fHistDataPM->GetMaximum();
+  Float_t minX = 1.001*fHistDataPM->GetXaxis()->GetXmin();
+  Float_t maxX = 0.999*fHistDataPM->GetXaxis()->GetXmax();
+  Int_t binSize = Int_t(1000*fHistDataPM->GetBinWidth(1));   // in MeV
+  Float_t minMinY = fHistSignal->GetMinimum();
+
   TCanvas *cSub = new TCanvas("cSub","signal, background subtracted",1400,1000);
-  cSub->Divide(2,2);
+  cSub->SetLeftMargin(0.15);
+  cSub->SetRightMargin(0.0);
+  cSub->SetTopMargin(0.002);
+  cSub->SetBottomMargin(0.0);
+  cSub->Divide(2,2,0.,0.);
   cSub->Draw();
 
-  Double_t minX  = fSignPM->GetBinLowEdge(1);       // minimum X value in axis
-  Double_t maxX  = fSignPM->GetBinLowEdge(fBins+1); // maximum X value in axis
-  if ( TMath::Abs(fDrawMin)<1.e-30 ) fDrawMin = minX;
-  if ( TMath::Abs(fDrawMax)<1.e-30 ) fDrawMax = maxX;
+  TVirtualPad* pad = cSub->cd(1);
+  pad->SetLeftMargin(0.15);
+  pad->SetRightMargin(0.0);
+  pad->SetTopMargin(0.005);
+  pad->SetBottomMargin(0.0);
+  TH2F *range1=new TH2F("range1","",10,minX,maxX,10,minY,maxY);
+  range1->SetStats(kFALSE);
+  range1->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
+  range1->GetYaxis()->CenterTitle();
+  range1->GetYaxis()->SetLabelSize(0.05);
+  range1->GetYaxis()->SetTitleSize(0.06);
+  range1->GetYaxis()->SetTitleOffset(0.8);
+  range1->Draw();
+  fHistDataPM->SetLineColor(1);
+  fHistDataPM->SetLineWidth(2);
+  //  fHistDataPM->SetMarkerStyle(21);
+  fHistDataPM->Draw("Psame");
+  TLatex *latex = new TLatex();
+  latex->SetNDC();
+  latex->SetTextSize(0.05);
+  latex->DrawLatex(0.2, 0.95, "Background un-substracted");
+  TLine line;
+  line.SetLineWidth(1);
+  line.SetLineStyle(2);
+  line.DrawLine(fIntMin, minY, fIntMin, maxY);
+  line.DrawLine(fIntMax, minY, fIntMax, maxY);
+
+  pad = cSub->cd(2);
+  pad->SetLeftMargin(0.);
+  pad->SetRightMargin(0.005);
+  pad->SetTopMargin(0.005);
+  pad->SetBottomMargin(0.0);
+  TH2F *range2=new TH2F("range2","",10,minX,maxX,10,minY,maxY);
+  range2->SetStats(kFALSE);
+  range2->Draw();
+  fHistBackground->SetLineColor(4);
+  fHistBackground->SetLineWidth(2);
+  //  fHistBackground->SetMarkerColor(4);
+  //  fHistBackground->SetMarkerStyle(6);
+  fHistBackground->Draw("Psame");
+  latex->DrawLatex(0.05, 0.95, "Like-sign background");
+  line.DrawLine(fIntMin, minY, fIntMin, maxY);
+  line.DrawLine(fIntMax, minY, fIntMax, maxY);
+  TLegend *legend = new TLegend(0.65, 0.70, 0.98, 0.98);
+  legend->SetFillColor(0);
+  legend->SetMargin(0.15);
+  legend->AddEntry(fHistDataPM, "N_{+-}", "l");
+  legend->AddEntry(fHistDataPP, "N_{++}", "l");
+  legend->AddEntry(fHistDataMM, "N_{--}", "l");
+  legend->AddEntry(fHistSignal, "N_{+-} - 2 #sqrt{N_{++} #times N_{--}}", "l");
+  legend->AddEntry(fHistBackground, "2 #sqrt{N_{++} #times N_{--}}", "l");
+  legend->Draw();
+
   
-  cSub->cd(4);
-  fSignPM->GetXaxis()->SetRangeUser(fDrawMin, fDrawMax);
-  fSignPM->SetLineColor(1);
-  fSignPM->SetLineWidth(2);
-  fSignPM->SetMarkerStyle(6);
-  fSignPM->DrawCopy("P");
-
-  fBackground->SetLineColor(4);
-  fBackground->SetMarkerColor(4);
-  fBackground->SetMarkerStyle(6);
-  fBackground->DrawCopy("Psame");
-
-  fSignal->SetMarkerStyle(20);
-  fSignal->SetMarkerColor(2);
-  fSignal->DrawCopy("Psame");
-
-  cSub->cd(1);
-  fSignal->DrawCopy("P");
-
-  cSub->cd(2);
-  fSignPM->DrawCopy("P");
-
-  cSub->cd(3);
-  fSignPP->DrawCopy("P");
+  pad = cSub->cd(3);
+  pad->SetLeftMargin(0.15);
+  pad->SetRightMargin(0.0);
+  pad->SetTopMargin(0.0);
+  pad->SetBottomMargin(0.15);
+  TH2F *range3=new TH2F("range3","",10,minX,maxX,10,minMinY,maxY);
+  range3->SetStats(kFALSE);
+  range3->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
+  range3->GetYaxis()->CenterTitle();
+  range3->GetYaxis()->SetLabelSize(0.05);
+  range3->GetYaxis()->SetTitleSize(0.06);
+  range3->GetYaxis()->SetTitleOffset(0.8);
+  range3->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
+  range3->GetXaxis()->CenterTitle();
+  range3->GetXaxis()->SetLabelSize(0.05);
+  range3->GetXaxis()->SetTitleSize(0.06);
+  range3->GetXaxis()->SetTitleOffset(1.0);
+  range3->Draw();
+  fHistDataPM->Draw("Psame");
+  fHistDataPP->SetLineWidth(2);
+  fHistDataPP->SetLineColor(6);
+  fHistDataMM->SetLineWidth(2);
+  fHistDataMM->SetLineColor(8);
+  fHistDataPP->Draw("Psame");
+  fHistDataMM->Draw("Psame");
+  line.DrawLine(minX, 0.,maxX, 0.);
+  line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
+  line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
+
+  pad = cSub->cd(4);
+  pad->SetLeftMargin(0.0);
+  pad->SetRightMargin(0.005);
+  pad->SetTopMargin(0.0);
+  pad->SetBottomMargin(0.15);
+  TH2F *range4=new TH2F("range4","",10,minX,maxX,10,minMinY,maxY);
+  range4->SetStats(kFALSE);
+  range4->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
+  range4->GetXaxis()->CenterTitle();
+  range4->GetXaxis()->SetLabelSize(0.05);
+  range4->GetXaxis()->SetTitleSize(0.06);
+  range4->GetXaxis()->SetTitleOffset(1.0);
+  range4->Draw();
+  fHistSignal->SetLineWidth(2);
+  fHistSignal->SetLineColor(2);
+  fHistSignal->Draw("Psame");
+  latex->DrawLatex(0.05, 0.95, "Like-sign background substracted");
+  if(fProcessed) DrawStats(0.05, 0.6, 0.5, 0.9);
+  line.DrawLine(minX, 0.,maxX, 0.);
+  line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
+  line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
 }
 
index ea36afa..06f836a 100644 (file)
 //#                                                           #
 //#############################################################
 
+/*
+  Class used for extracting the signal from an invariant mass spectrum.
+  It implements the AliDielectronSignalBase class and uses the like-sign
+  substraction method for estimating the signal and background.
+  There is no fitting in this class, only bin counting.
+
+  Example usage:
+   AliDielectronSignalExt *signalProcess = new AliDielectronSignalExt();
+   TObjArray *histoArray = new TObjArray();
+   histoArray->Add(signalPP);                  // the order of putting the histograms in the array is important!!
+   histoArray->Add(signalPM);
+   histoArray->Add(signalMM);
+   signalProcess->SetMethod(AliDielectronSignalBase::kLikeSign);  // or kEventMixing
+   signalProcess->SetIntegralRange(3.0,3.15);   // J/Psi peak
+   signalProcess->SetRebin(2);                  // rebin the histograms
+   signalProcess->Process(histoArray);
+   signalProcess->Draw("stat");
+   signalProcess->Print();
+
+*/
+
 #include <TVectorT.h>
 #include <TString.h>
 #include <TH1.h>
@@ -28,6 +49,7 @@
 class AliDielectronSignalExt : public AliDielectronSignalBase {
 
 public:
   AliDielectronSignalExt();
   AliDielectronSignalExt(const char*name, const char* title);
 
@@ -37,46 +59,14 @@ public:
   void ProcessLS(TObjArray* const arrhist);  // like-sign method
   void ProcessEM(TObjArray* const arrhist);  // event mixing method
   
-  // getters
-  TH1F* GetHistogramSignal()     const { return fSignal; } 
-  TH1F* GetHistogramBackground() const { return fBackground; }
-
-  // setters
-  void SetMethod(Int_t method){ fMethod=method; }
-  void SetHistograms(TH1F* const unlike, TH1F* const backg, TH1F* const signal);
-  void SetRebin(Int_t rebin) {fRebin = rebin;}
-  void SetDrawRange(Double_t min, Double_t max){fDrawMin=min; fDrawMax=max;}
-  void SetFitRange(Double_t min, Double_t max) {fFitMin=min;fFitMax=max;}
-
-  Double_t GetFitMin()      const { return fFitMin; }
-  Double_t GetFitMax()      const { return fFitMax; }
-  void Rebin(Int_t rebin);
-  
   virtual void Draw(const Option_t* option = "");
 
-
 private:
-  
-  TH1F *fSignPM;              // histogram of unlike sign (plus-minus)
-  TH1F *fSignPP;              // histogram of like sign (plus-plus)
-  TH1F *fSignMM;              // histogram of like sign (minus-minus)
-  TH1F *fBackground;          // histogram of like-sign background
-  TH1F *fSignal;              // histogram of subtracted signal
-  
-  Int_t    fMethod;           // subtraction method. 1(like-sign), 2(event mixing)
-  Int_t    fRebin;            // number of histogram rebin iteration
-  Int_t    fBins;             // number of bins in X axis
-  Double_t fDrawMin;          // minimum X when drawing 
-  Double_t fDrawMax;          // maximum X when drawing
-  Double_t fFitMin;           // fit range min
-  Double_t fFitMax;           // fit range max
 
   AliDielectronSignalExt(const AliDielectronSignalExt &c);
   AliDielectronSignalExt &operator=(const AliDielectronSignalExt &c);
 
-  ClassDef(AliDielectronSignalExt,1)         // Dielectron SignalFunc
+  ClassDef(AliDielectronSignalExt,2)    // Dielectron SignalFunc
 };
 
-
-
 #endif
index eeef9d3..11069cb 100644 (file)
@@ -34,6 +34,9 @@ can be used.
 #include <TMath.h>
 #include <TString.h>
 #include <TPaveText.h>
+#include <TList.h>
+#include <TFitResult.h>
+#include <../hist/hist/src/TF1Helper.h>
 
 #include <AliLog.h>
 
@@ -43,17 +46,13 @@ ClassImp(AliDielectronSignalFunc)
 
 AliDielectronSignalFunc::AliDielectronSignalFunc() :
   AliDielectronSignalBase(),
-  fSignal(0x0),
-  fBackground(0x0),
-  fSigBack(0x0),
-  fVInitParam(0),
-  fFitOpt("MNQE"),
-  fUseIntegral(kFALSE),
-  fFitMin(2.5),
-  fFitMax(4),
-  fParM(-1),
-  fParMres(-1),
-  fSignalHist(0x0)
+  fFuncSignal(0x0),
+  fFuncBackground(0x0),
+  fFuncSigBack(0x0),
+  fParMass(1),
+  fParMassWidth(2),
+  fFitOpt("SMNQE"),
+  fUseIntegral(kFALSE)
 {
   //
   // Default Constructor
@@ -64,17 +63,13 @@ AliDielectronSignalFunc::AliDielectronSignalFunc() :
 //______________________________________________
 AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) :
   AliDielectronSignalBase(name, title),
-  fSignal(0x0),
-  fBackground(0x0),
-  fSigBack(0x0),
-  fVInitParam(0),
-  fFitOpt("MNQ"),
-  fUseIntegral(kFALSE),
-  fFitMin(2.5),
-  fFitMax(4),
-  fParM(-1),
-  fParMres(-1),
-  fSignalHist(0x0)
+  fFuncSignal(0x0),
+  fFuncBackground(0x0),
+  fFuncSigBack(0x0),
+  fParMass(1),
+  fParMassWidth(2),
+  fFitOpt("SMNQE"),
+  fUseIntegral(kFALSE)
 {
   //
   // Named Constructor
@@ -87,98 +82,248 @@ AliDielectronSignalFunc::~AliDielectronSignalFunc()
   //
   // Default Destructor
   //
-  if (fSigBack) delete fSigBack;
+  if(fFuncSignal) delete fFuncSignal;
+  if(fFuncBackground) delete fFuncBackground;
+  if(fFuncSigBack) delete fFuncSigBack;
 }
 
 
 //______________________________________________
-void AliDielectronSignalFunc::Process(TH1 * const hist)
+void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
 {
   //
-  // Fit the mass histogram with the function and retrieve the parameters
+  // Fit the invariant mass histograms and retrieve the signal and background
   //
-  
-  //reset result arrays
-  Reset();
+  switch(fMethod) {
+  case kFitted :
+    ProcessFit(arrhist);
+    break;
 
-  //sanity check
-  if (!fSigBack) {
-    AliError("Use 'SetFunctions(TF1*,TF1*)' or 'SetDefaults(Int_t)' to setup the signal functions first'!");
-    return;
+  case kLikeSign :
+    ProcessLS(arrhist);
+    break;
+
+  case kEventMixing :
+      //ProcessEM(arrhist); //yet to be implemented...
+    break;
+
+  default :
+    AliError("Background substraction method not known!");
   }
-  
-  //set current signal hist pointer
-  fSignalHist=hist;
-  
-  //initial the fit function to its default parameters
-  if (fVInitParam.GetMatrixArray()) fSigBack->SetParameters(fVInitParam.GetMatrixArray());
+}
 
-  //Perform fit and retrieve values
-  hist->Fit(fSigBack,fFitOpt.Data(),"",fFitMin,fFitMax);
 
-  //retrieve values and errors
-  //TODO: perhpas implement different methods to retrieve the valus
-  fSignal->SetParameters(fSigBack->GetParameters());
-  fBackground->SetParameters(fSigBack->GetParameters()+fSignal->GetNpar());
+//______________________________________________
+void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
+  //
+  // Fit the +- invariant mass distribution only
+  // Here we assume that the combined fit function is a sum of the signal and background functions
+  //    and that the signal function is always the first term of this sum
+  //
 
-  //TODO: proper error estimate?!?
-  //      perhaps this is not possible in a general way?
+  fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
+  fHistDataPM->Sumw2();
+  if(fRebin>1)
+    fHistDataPM->Rebin(fRebin);
+
+  fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal", 
+                        fHistDataPM->GetXaxis()->GetNbins(),
+                        fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
+  fHistBackground = new TH1F("HistBackground", "Like-sign contribution", 
+                            fHistDataPM->GetXaxis()->GetNbins(),
+                            fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
   
-  // signal and background histograms
-  TH1 *hSignal=(TH1*)hist->Clone("DieleSignalHist");
-  TH1 *hBackground=(TH1*)hist->Clone("DieleBackgroundHist");
-
-  fSignal->SetRange(hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax());
-  fBackground->SetRange(hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax());
+  // the starting parameters of the fit function and their limits can be tuned 
+  // by the user in its macro
+  fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
+  TFitResultPtr pmFitPtr = fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
+  TFitResult *pmFitResult = pmFitPtr.Get();  
+  fFuncSignal->SetParameters(fFuncSigBack->GetParameters());
+  fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar());
   
-  hSignal->Add(fBackground,-1);
-  hBackground->Add(fSignal,-1);
+  for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
+    Double_t m = fHistDataPM->GetBinCenter(iBin);
+    Double_t pm = fHistDataPM->GetBinContent(iBin);
+    Double_t epm = fHistDataPM->GetBinError(iBin);
+    Double_t bknd = fFuncBackground->Eval(m);
+    Double_t ebknd = 0;
+    for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
+      for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
+        TF1 gradientIpar("gradientIpar",
+       ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
+        TF1 gradientJpar("gradientJpar",
+       ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
+       ebknd += pmFitResult->CovMatrix(iPar,jPar)*
+                gradientIpar.Eval(m)*gradientJpar.Eval(m)*
+                (iPar==jPar ? 1.0 : 2.0);
+      }
+    }
+    Double_t signal = pm-bknd;
+    Double_t error = TMath::Sqrt(epm*epm+ebknd);
+    fHistSignal->SetBinContent(iBin, signal);
+    fHistSignal->SetBinError(iBin, error);
+    fHistBackground->SetBinContent(iBin, bknd);
+    fHistBackground->SetBinError(iBin, TMath::Sqrt(ebknd));
+  }
 
-  
-  //get values and errors signal and background
-  Double_t signal=0, signal_err=0;
-  Double_t background=0, background_err=0;
-  Double_t sigPlusBack=0, sigPlusBack_err=0;
-
-  if (!fUseIntegral){
-    signal=hSignal->IntegralAndError(hSignal->FindBin(GetIntegralMin()),
-                                     hSignal->FindBin(GetIntegralMax()), signal_err);
-    background=hBackground->IntegralAndError(hBackground->FindBin(GetIntegralMin()),
-                                             hBackground->FindBin(GetIntegralMax()), background_err);
-    sigPlusBack=hist->IntegralAndError(hist->FindBin(GetIntegralMin()),
-                                       hist->FindBin(GetIntegralMax()), sigPlusBack_err);
-  } else {
-    signal=fSignal->Integral(GetIntegralMin(),GetIntegralMax())/hist->GetBinWidth(hist->FindBin((GetIntegralMax()-GetIntegralMin())/2.));
-    background=fBackground->Integral(GetIntegralMin(),GetIntegralMax())/hist->GetBinWidth(hist->FindBin((GetIntegralMax()-GetIntegralMin())/2.));
-    sigPlusBack=fSigBack->Integral(GetIntegralMin(),GetIntegralMax())/hist->GetBinWidth(hist->FindBin((GetIntegralMax()-GetIntegralMin())/2.));
-    //TODO: Error calculation
+  if(fUseIntegral) {
+    // signal
+    fValues(0) = fFuncSignal->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
+    fErrors(0) = 0;
+    for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
+      for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) {
+        TF1 gradientIpar("gradientIpar",
+       ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0);
+        TF1 gradientJpar("gradientJpar",
+       ROOT::TF1Helper::TGradientParFunction(jPar,fFuncSignal),0,0,0);
+       fErrors(0) += pmFitResult->CovMatrix(iPar,jPar)*
+         gradientIpar.Integral(fIntMin,fIntMax)*gradientJpar.Integral(fIntMin,fIntMax)*
+         (iPar==jPar ? 1.0 : 2.0);
+      }
+    } 
+    // background
+    fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
+    fErrors(1) = 0;
+    for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
+      for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
+        TF1 gradientIpar("gradientIpar",
+       ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
+        TF1 gradientJpar("gradientJpar",
+       ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
+       fErrors(1) += pmFitResult->CovMatrix(iPar,jPar)*
+         gradientIpar.Integral(fIntMin, fIntMax)*gradientJpar.Integral(fIntMin, fIntMax)*
+         (iPar==jPar ? 1.0 : 2.0);
+      }
+    }
   }
-  
-  //set values
-  SetSignal(signal,signal_err);
-  SetBackground(background,background_err);
-  SetSignificanceAndSOB();
-  if (fParM>-1){
-    SetMass(fSigBack->GetParameter(fParM), fSigBack->GetParError(fParM));
-    SetMassWidth(fSigBack->GetParameter(fParMres), fSigBack->GetParError(fParMres));
+  else {
+    // signal
+    fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
+                                           fHistSignal->FindBin(fIntMax), fErrors(0));
+    // background
+    fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
+                                                   fHistBackground->FindBin(fIntMax), 
+                                                   fErrors(1));
   }
-  //cleanup
-  delete hSignal;
-  delete hBackground;
+  // S/B and significance
+  SetSignificanceAndSOB();
+  fValues(4) = fFuncSigBack->GetParameter(fParMass);
+  fErrors(4) = fFuncSigBack->GetParError(fParMass);
+  fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
+  fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
+
+  fProcessed = kTRUE;
 }
 
 //______________________________________________
-void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
-{
+void AliDielectronSignalFunc::ProcessLS(TObjArray * const arrhist) {
   //
-  // Fit the mass histogram with the function and retrieve the parameters
+  // Substract background using the like-sign spectrum
   //
-  
-  TH1 *hist=(TH1*)arrhist->At(1);
-  Process(hist);
+  fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP");  // ++    SE
+  fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
+  fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM");  // --    SE 
+  if (fRebin>1) { 
+    fHistDataPP->Rebin(fRebin);
+    fHistDataPM->Rebin(fRebin);
+    fHistDataMM->Rebin(fRebin);
+  }   
+  fHistDataPP->Sumw2();
+  fHistDataPM->Sumw2();
+  fHistDataMM->Sumw2();
+
+  fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal", 
+                        fHistDataPM->GetXaxis()->GetNbins(),
+                        fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
+  fHistBackground = new TH1F("HistBackground", "Like-sign contribution", 
+                            fHistDataPM->GetXaxis()->GetNbins(),
+                            fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
+
+  // fit the +- mass distribution
+  fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
+  fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
+  // declare the variables where the like-sign fit results will be stored
+  TFitResult *ppFitResult = 0x0;
+  TFitResult *mmFitResult = 0x0;
+  // fit the like sign background
+  TF1 *funcClonePP = (TF1*)fFuncBackground->Clone("funcClonePP");
+  TF1 *funcCloneMM = (TF1*)fFuncBackground->Clone("funcCloneMM");
+  fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
+  TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
+  ppFitResult = ppFitPtr.Get();
+  fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
+  TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
+  mmFitResult = mmFitPtr.Get();
+
+  for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
+    Double_t m = fHistDataPM->GetBinCenter(iBin);
+    Double_t pm = fHistDataPM->GetBinContent(iBin);
+    Double_t pp = funcClonePP->Eval(m);
+    Double_t mm = funcCloneMM->Eval(m);
+    Double_t epm = fHistDataPM->GetBinError(iBin);
+    Double_t epp = 0;
+    for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
+      for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
+        TF1 gradientIpar("gradientIpar",
+       ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
+        TF1 gradientJpar("gradientJpar",
+       ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
+       epp += ppFitResult->CovMatrix(iPar,jPar)*
+              gradientIpar.Eval(m)*gradientJpar.Eval(m)*
+              (iPar==jPar ? 1.0 : 2.0);
+      }
+    }
+    Double_t emm = 0;
+    for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
+      for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
+        TF1 gradientIpar("gradientIpar",
+       ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
+        TF1 gradientJpar("gradientJpar",
+       ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
+       emm += mmFitResult->CovMatrix(iPar,jPar)*
+              gradientIpar.Eval(m)*gradientJpar.Eval(m)*
+              (iPar==jPar ? 1.0 : 2.0);
+      }
+    }
+
+    Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
+    Double_t background = 2.0*TMath::Sqrt(pp*mm);
+    // error propagation on the signal calculation above
+    Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
+    Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
+    fHistSignal->SetBinContent(iBin, signal);
+    fHistSignal->SetBinError(iBin, esignal);
+    fHistBackground->SetBinContent(iBin, background);
+    fHistBackground->SetBinError(iBin, ebackground);
+  }
+
+  // signal
+  fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
+                                         fHistSignal->FindBin(fIntMax), fErrors(0));
+  // background
+  fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
+                                                 fHistBackground->FindBin(fIntMax), 
+                                                 fErrors(1));
+  // S/B and significance
+  SetSignificanceAndSOB();
+  fValues(4) = fFuncSigBack->GetParameter(fParMass);
+  fErrors(4) = fFuncSigBack->GetParError(fParMass);
+  fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
+  fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
+
+  fProcessed = kTRUE;
 }
+
+//______________________________________________
+//void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
+  //
+  // Substract background with the event mixing technique
+  //
+  //AliError("Event mixing for background substraction method not implemented!");
+//}
+
 //______________________________________________
-void AliDielectronSignalFunc::SetFunctions(TF1 * const sig, TF1 * const back, TF1 * const combined,
+void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
                                            Int_t parM, Int_t parMres)
 {
   //
@@ -186,30 +331,16 @@ void AliDielectronSignalFunc::SetFunctions(TF1 * const sig, TF1 * const back, TF
   // Note: The process method assumes that the first n parameters in the
   //       combined fit function correspond to the n parameters of the signal function
   //       and the n+1 to n+m parameters to the m parameters of the background function!!!
-  // if parM and parMres are set this information can be used in the drawing function
-  //   to also show in the statistics box the mass and mass resolution
-
+  
   if (!sig||!back||!combined) {
     AliError("Both, signal and background function need to be set!");
     return;
   }
-  fSignal=sig;
-  fBackground=back;
-  fSigBack=combined;
-
-  InitParams();
-  fParM=parM;
-  fParMres=parMres;
-}
-
-//______________________________________________
-void AliDielectronSignalFunc::InitParams()
-{
-  //
-  // initilise common fit parameters
-  //
-  fVInitParam.ResizeTo(fSigBack->GetNpar());
-  fVInitParam.SetElements(fSigBack->GetParameters());
+  fFuncSignal=sig;
+  fFuncBackground=back;
+  fFuncSigBack=combined;
+  fParMass=parM;
+  fParMassWidth=parMres;
 }
 
 //______________________________________________
@@ -223,49 +354,41 @@ void AliDielectronSignalFunc::SetDefaults(Int_t type)
   // type = 3: Crystal-Ball function
   // type = 4: Crystal-Ball signal + exponential background
   //
-
   
   if (type==0){
-    fSignal=new TF1("DieleSignal","gaus",2.5,4);
-    fBackground=new TF1("DieleBackground","pol1",2.5,4);
-    fSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
-    
-    fSigBack->SetParameters(1,3.1,.05,2.5,1);
-    fSigBack->SetParLimits(0,0,10000000);
-    fSigBack->SetParLimits(1,3.05,3.15);
-    fSigBack->SetParLimits(2,.02,.1);
-    
-    SetFunctions(fSignal,fBackground,fSigBack,1,2);
-
-  } else if (type==1){
-    fSignal=new TF1("DieleSignal","gaus",2.5,4);
-    fBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
-    fSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
+    fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
+    fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
+    fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
     
-    fSigBack->SetParameters(1,3.1,.05,1,2.5,1);
-    fSigBack->SetParLimits(0,0,10000000);
-    fSigBack->SetParLimits(1,3.05,3.15);
-    fSigBack->SetParLimits(2,.02,.1);
+    fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
+    fFuncSigBack->SetParLimits(0,0,10000000);
+    fFuncSigBack->SetParLimits(1,3.05,3.15);
+    fFuncSigBack->SetParLimits(2,.02,.1);
+  } 
+  else if (type==1){
+    fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
+    fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
+    fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
     
-    SetFunctions(fSignal,fBackground,fSigBack,1,2);
-
-  } else if (type==2){
+    fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
+    fFuncSigBack->SetParLimits(0,0,10000000);
+    fFuncSigBack->SetParLimits(1,3.05,3.15);
+    fFuncSigBack->SetParLimits(2,.02,.1);
+  } 
+  else if (type==2){
     // half gaussian, half exponential signal function
     // exponential background
-    fSignal = new     TF1("DieleSignal","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",2.5,4);
-    fBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
-    fSigBack=new TF1("DieleCombined","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))+[4]*exp(-(x-[5])/[6])+[7]",2.5,4);
-    fSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
+    fFuncSignal = new TF1("DieleSignal","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",2.5,4);
+    fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
+    fFuncSigBack = new TF1("DieleCombined","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))+[4]*exp(-(x-[5])/[6])+[7]",2.5,4);
+    fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
     
-    fSigBack->SetParLimits(0,0,10000000);
-    fSigBack->SetParLimits(1,3.05,3.15);
-    fSigBack->SetParLimits(2,.02,.1);
-    fSigBack->FixParameter(6,2.5);
-    fSigBack->FixParameter(7,0);
-    SetFunctions(fSignal,fBackground,fSigBack,1,2);
-
+    fFuncSigBack->SetParLimits(0,0,10000000);
+    fFuncSigBack->SetParLimits(1,3.05,3.15);
+    fFuncSigBack->SetParLimits(2,.02,.1);
+    fFuncSigBack->FixParameter(6,2.5);
+    fFuncSigBack->FixParameter(7,0);
   }
-  
 }
 
 
@@ -275,22 +398,22 @@ void AliDielectronSignalFunc::Draw(const Option_t* option)
   //
   // Draw the fitted function
   //
-
+  
   TString drawOpt(option);
   drawOpt.ToLower();
 
   Bool_t optStat=drawOpt.Contains("stat");
   
-  fSigBack->SetNpx(200);
-  fSigBack->SetRange(GetIntegralMin(),GetIntegralMax());
-  fBackground->SetNpx(200);
-  fBackground->SetRange(GetIntegralMin(),GetIntegralMax());
+  fFuncSigBack->SetNpx(200);
+  fFuncSigBack->SetRange(fIntMin,fIntMax);
+  fFuncBackground->SetNpx(200);
+  fFuncBackground->SetRange(fIntMin,fIntMax);
   
-  TGraph *grSig=new TGraph(fSigBack);
+  TGraph *grSig=new TGraph(fFuncSigBack);
   grSig->SetFillColor(kGreen);
   grSig->SetFillStyle(3001);
 
-  TGraph *grBack=new TGraph(fBackground);
+  TGraph *grBack=new TGraph(fFuncBackground);
   grBack->SetFillColor(kRed);
   grBack->SetFillStyle(3001);
 
@@ -300,12 +423,12 @@ void AliDielectronSignalFunc::Draw(const Option_t* option)
   grBack->SetPoint(0,grBack->GetX()[0],0.);
   grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
   
-  fSigBack->SetRange(fFitMin,fFitMax);
-  fBackground->SetRange(fFitMin,fFitMax);
+  fFuncSigBack->SetRange(fFitMin,fFitMax);
+  fFuncBackground->SetRange(fFitMin,fFitMax);
   
   if (!drawOpt.Contains("same")){
-    if (fSignalHist){
-      fSignalHist->Draw();
+    if (fHistDataPM){
+      fHistDataPM->Draw();
       grSig->Draw("f");
     } else {
       grSig->Draw("af");
@@ -313,11 +436,21 @@ void AliDielectronSignalFunc::Draw(const Option_t* option)
   } else {
     grSig->Draw("f");
   }
-  grBack->Draw("f");
-  fSigBack->Draw("same");
-  fBackground->Draw("same");
+  if(fMethod==kFitted) grBack->Draw("f");
+  fFuncSigBack->Draw("same");
+  fFuncSigBack->SetLineWidth(2);
+  if(fMethod==kLikeSign) {
+    fHistDataPP->SetLineWidth(2);
+    fHistDataPP->SetLineColor(6);
+    fHistDataPP->Draw("same");
+    fHistDataMM->SetLineWidth(2);
+    fHistDataMM->SetLineColor(8); 
+    fHistDataMM->Draw("same");
+  }
+
+  if(fMethod==kFitted) 
+    fFuncBackground->Draw("same");
 
   if (optStat) DrawStats();
-    
+  
 }
-
index b3da553..22b8690 100644 (file)
 //#                                                           #
 //#############################################################
 
+/*
+  Class used for extracting the signal from an invariant mass spectrum.
+  It implements the AliDielectronSignalBase class and it uses user provided
+  functions to fit the unlike-sign spectrum (and the like-sign one).
+  
+  Example usage:
+
+  AliDielectronSignalFunc *signalProcess = new AliDielectronSignalFunc();
+  TObjArray *histoArray = new TObjArray();
+  histoArray->Add(signalPP);            // add the spectrum histograms to the array
+  histoArray->Add(signalPM);            // the order is important !!!
+  histoArray->Add(signalMM);
+  // set the extraction method 
+  // AliDielectronSignalBase::kFitted       -->  fit only the unlike-sign spectrum and extract everything from that
+  // AliDielectronSignalBase::kLikeSign     -->  fit both the unlike- and like-sign spectra
+  // AliDielectronSignalBase::kEventMixing  -->  fit both the unlike- and like-sign spectra from event mixing
+  signalProcess->SetMethod(AliDielectronSignalBase::kLikeSign); 
+  // Initialize the functions to be used and pass them to the signal object
+  // External preparation of the functions can(should) be done as this can be a 5 or more parameter fit
+  TF1* gaus = new TF1("gaus", "gaus", 0., 4.);
+  TF1* expo = new TF1("expo", "[0]*exp([1]*x)", 0., 4.);
+  TF1* combined = new TF1("combined", "gaus + [3]*exp([4]*x)", 0.,4.);
+  combined->SetParameter(1, 3.1);
+  combined->SetParameter(1, 0.1);
+  signalPP->Fit(expo, "SME", "", 2.4, 4.0);
+  signalPM->Fit(gaus, "SME", "", 3.0, 3.15);
+  Double_t pars[5];
+  gaus->GetParameters(&pars[0]);
+  expo->GetParameters(&pars[3]);
+  combined->SetParameters(pars);
+  combined->SetParLimits(1, 3.05, 3.15);
+  combined->SetParLimits(2, 0.03, 0.1);
+  signalProcess->SetFunctions(combined, gaus, expo, 1, 2);
+
+  signalProcess->SetFitRange(2.4,4.0);
+  // Use the integral of the fit function to estimate the signal or not
+  // The background will always be estimated from the fit
+  //  signalProcess->SetUseIntegral(kTRUE);  
+  signalProcess->SetFitOption("SME");
+  // Give the range where the signal is calculated
+  signalProcess->SetIntegralRange(3.0,3.15);
+  signalProcess->SetRebin(2);
+  signalProcess->Process(histoArray);
+  signalProcess->Draw("stat");
+  signalProcess->Print();
+*/
+
+
 #include <TVectorT.h>
 #include <TString.h>
+#include <TH1F.h>
 
 #include "AliDielectronSignalBase.h"
 
@@ -28,52 +77,43 @@ class AliDielectronSignalFunc : public AliDielectronSignalBase {
 public:
   AliDielectronSignalFunc();
   AliDielectronSignalFunc(const char*name, const char* title);
+  AliDielectronSignalFunc(const AliDielectronSignalFunc &c);
+  AliDielectronSignalFunc &operator=(const AliDielectronSignalFunc &c);
 
   virtual ~AliDielectronSignalFunc();
 
-  void Process(TH1 * const hist);
   virtual void Process(TObjArray * const arrhist);
+  void ProcessFit(TObjArray * const arrhist);      // fit the SE +- distribution
+  void ProcessLS(TObjArray * const arrhist);       // substract the fitted SE like-sign background
+  //void ProcessEM(TObjArray * const arrhist);       // substract the fitted SE+ME like-sign background ...not yet implemented
   
-  void SetFunctions(TF1 * const sig, TF1 * const back, TF1 * const combined, Int_t parM=-1, Int_t parMres=-1);
-  void InitParams();
-  void SetFitOption(const char* opt) {fFitOpt=opt;}
+  void SetUseIntegral(Bool_t flag=kTRUE) {fUseIntegral = flag;};
+  void SetFunctions(TF1 * const combined, TF1 * const sig=0, TF1 * const back=0, Int_t parM=1, Int_t parMres=2);
+  void SetFitOption(const char* opt) {
+    fFitOpt=opt; 
+    fFitOpt.ToLower(); 
+    if(!fFitOpt.Contains("s")) fFitOpt += "s";
+  }
   void SetDefaults(Int_t type);
-  void SetFitRange(Double_t min, Double_t max) {fFitMin=min;fFitMax=max;}
-
-  void SetUseIntegral(Bool_t integral=kTRUE) { fUseIntegral=integral; }
-  
-  const char* GetFitOption()          const { return fFitOpt.Data(); }
-  TF1*  GetSignalFunction()     const { return fSignal;        }
-  TF1*  GetBackgroundFunction() const { return fBackground;    }
-  TF1*  GetCombinedFunction()   const { return fSigBack;       }
-  
-  Double_t GetFitMin()      const { return fFitMin; }
-  Double_t GetFitMax()      const { return fFitMax; }
+    
+  TF1*  GetSignalFunction()     { return fFuncSignal;        }
+  TF1*  GetBackgroundFunction() { return fFuncBackground;    }
+  TF1*  GetCombinedFunction()   { return fFuncSigBack;       }
   
   virtual void Draw(const Option_t* option = "");
   
 private:
-  TF1 *fSignal;                // Function for the signal description
-  TF1 *fBackground;            // Function for the background description
-  TF1 *fSigBack;               // Combined function signal plus background
-  TVectorD fVInitParam;        // initial fit parameters
-  TString fFitOpt;             // fit option used
-
-  Bool_t fUseIntegral;         // use integral instead of bin counts
-  
-  Double_t fFitMin;            // fit range min
-  Double_t fFitMax;            // fit range max
-
-  Int_t fParM;                 // Paramter which defines the mass
-  Int_t fParMres;              // Paramter which defines the resolution of the mass
 
-  TH1 *fSignalHist;            // Current signal histogram
+  TF1 *fFuncSignal;                // Function for the signal description
+  TF1 *fFuncBackground;            // Function for the background description
+  TF1 *fFuncSigBack;               // Combined function signal plus background
+  Int_t fParMass;                  // the index of the parameter corresponding to the resonance mass
+  Int_t fParMassWidth;             // the index of the parameter corresponding to the resonance mass width
   
-  AliDielectronSignalFunc(const AliDielectronSignalFunc &c);
-  AliDielectronSignalFunc &operator=(const AliDielectronSignalFunc &c);
+  TString fFitOpt;             // fit option used
+  Bool_t fUseIntegral;         // use the integral of the fitted functions to extract signal and background
 
-  ClassDef(AliDielectronSignalFunc,1)         // Dielectron SignalFunc
+  ClassDef(AliDielectronSignalFunc,2)         // Dielectron SignalFunc
 };
 
-
 #endif
index 31a73b4..ed4ae10 100644 (file)
@@ -226,8 +226,8 @@ void AliDielectronSpectrum::CreateCFSpectrum()
   for (Int_t iMethod=0; iMethod<fSignalMethods.GetEntries(); ++iMethod){
     TString name(fSignalMethods.UncheckedAt(iMethod)->GetName());
     if (fStepSignal){
-      fCFSpectrum->SetStepTitle(steps++,(name+" (Siganl)").Data());
-      fCFSpectrum->SetStepTitle(steps++,(name+" (Corrected Siganl)").Data());
+      fCFSpectrum->SetStepTitle(steps++,(name+" (Signal)").Data());
+      fCFSpectrum->SetStepTitle(steps++,(name+" (Corrected Signal)").Data());
     }
     if (fStepSignificance){
       fCFSpectrum->SetStepTitle(steps++,(name+" (Significance)").Data());
index 12f521a..dec472c 100644 (file)
@@ -7,7 +7,7 @@ void InitCF();
 
 AliESDtrackCuts *SetupESDtrackCuts();
 
-TString names=("basicQ+SPDfirst+pt>.6+PID");
+TString names=("basicQ+SPDfirst+pt>.6+PID;basicQ+SPDany+pt>.6+PID");
 
 TObjArray *arrNames=names.Tokenize(";");
 
@@ -17,13 +17,13 @@ AliDielectron *fDiele=0x0;
 Int_t          fCutDefinition=0;
 Bool_t         fIsAOD=kFALSE;
 
-AliDielectron* ConfigJpsi2ee(Int_t fCutDefinition, Bool_t isAOD)
+AliDielectron* ConfigJpsi2ee(Int_t cutDefinition, Bool_t isAOD)
 {
   //
   // Setup the instance of AliDielectron
   //
   
-  fCutDefinition=fCutDefinition;
+  fCutDefinition=cutDefinition;
   fIsAOD=isAOD;
   
   // create the actual framework object
@@ -64,34 +64,35 @@ void SetupTrackCuts()
     fDiele->GetTrackFilter().AddCuts(SetupESDtrackCuts());
   } else {
     AliDielectronTrackCuts *trackCuts=new AliDielectronTrackCuts("trackCuts","trackCuts");
-    trackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
+    if (fCutDefinition==0)
+      trackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
+    else if (fCutDefinition==1)
+      trackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
     trackCuts->SetRequireTPCRefit(kTRUE);
     trackCuts->SetRequireITSRefit(kTRUE);
     fDiele->GetTrackFilter().AddCuts(trackCuts);
   }
 
-  if(fCutDefinition==0) {
-    //Pt cut ----------------------------------------------------------
-    AliDielectronVarCuts *pt = new AliDielectronVarCuts("ptCut","pt cut");
-    pt->AddCut(AliDielectronVarManager::kPt,0.6,1e30);
-
-    //AOD additions since there are no AliESDtrackCuts -----------------
-    //
-    if (fIsAOD){
-      // TPC #clusteres cut
-      pt->AddCut(AliDielectronVarManager::kNclsTPC,90.,160.);
-      pt->AddCut(AliDielectronVarManager::kEta,-0.8,0.8);
-      //TODO: DCA cuts to be investigated!!!
+  //Pt cut ----------------------------------------------------------
+  AliDielectronVarCuts *pt = new AliDielectronVarCuts("ptCut","pt cut");
+  pt->AddCut(AliDielectronVarManager::kPt,0.6,1e30);
+
+  //AOD additions since there are no AliESDtrackCuts -----------------
+  //
+  if (fIsAOD){
+    // TPC #clusteres cut
+    pt->AddCut(AliDielectronVarManager::kNclsTPC,90.,160.);
+    pt->AddCut(AliDielectronVarManager::kEta,-0.8,0.8);
+    //TODO: DCA cuts to be investigated!!!
 //       pt->AddCut(AliDielectronVarManager::kImpactParXY,-1.,1.);
 //       pt->AddCut(AliDielectronVarManager::kImpactParZ,-3.,3.);
-    }
-    fDiele->GetTrackFilter().AddCuts(pt);
-    
-    // PID cuts --------------------------------------------------------
-    AliDielectronPID *pid = new AliDielectronPID("PID10","TPC nSigma |e|<3 + |Pi|>3 + |P|>3 + TOF nSigma |e|<3");
-    pid->SetDefaults(10);
-    fDiele->GetTrackFilter().AddCuts(pid);
   }
+  fDiele->GetTrackFilter().AddCuts(pt);
+    
+  // PID cuts --------------------------------------------------------
+  AliDielectronPID *pid = new AliDielectronPID("PID10","TPC nSigma |e|<3 + |Pi|>3 + |P|>3 + TOF nSigma |e|<3");
+  pid->SetDefaults(10);
+  fDiele->GetTrackFilter().AddCuts(pid);
 }
 
 //______________________________________________________________________________________
@@ -129,7 +130,10 @@ AliESDtrackCuts *SetupESDtrackCuts()
   esdTrackCuts->SetMinNClustersTPC(90);
   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
   
-  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
+  if (fCutDefinition==0)
+    esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
+  else if (fCutDefinition==1)
+    esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
   
   return esdTrackCuts;
 }
@@ -172,12 +176,12 @@ void InitHistograms()
   histos->UserHistogram("Track","dXY","dXY;dXY [cm];#tracks",500,-1.,1.,AliDielectronVarManager::kImpactParXY);
   histos->UserHistogram("Track","dZ","dZ;dZ [cm];#tracks",600,-3.,3.,AliDielectronVarManager::kImpactParZ);
   histos->UserHistogram("Track","Eta_Phi","Eta Phi Map; Eta; Phi;#tracks",
-                        200,-1,1,200,0,6.285,AliDielectronVarManager::kEta,AliDielectronVarManager::kPhi);
+                        100,-1,1,144,0,6.285,AliDielectronVarManager::kEta,AliDielectronVarManager::kPhi);
 
   histos->UserHistogram("Track","dEdx_P","dEdx;P [GeV];TPC signal (arb units);#tracks",
-                        400,0.2,20.,200,0.,200.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCsignal,kTRUE);
+                        200,0.2,20.,100,0.,200.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCsignal,kTRUE);
   histos->UserHistogram("Track","TPCnSigmaEle_P","TPC number of sigmas Electrons;P [GeV];TPC number of sigmas Electrons;#tracks",
-                        400,0.2,20.,200,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle,kTRUE);
+                        200,0.2,20.,100,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle,kTRUE);
       
   //add histograms to Pair classes
   histos->UserHistogram("Pair","InvMass","Inv.Mass;Inv. Mass [GeV];#pairs",
index 85e518c..73de2e9 100644 (file)
@@ -45,7 +45,7 @@ void SetupTrackCuts()
     fDiele->GetTrackFilter().AddCuts(SetupESDtrackCuts());
   } else {
     AliDielectronTrackCuts *trackCuts=new AliDielectronTrackCuts("trackCuts","trackCuts");
-    trackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
+    trackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
     trackCuts->SetRequireTPCRefit(kTRUE);
     trackCuts->SetRequireITSRefit(kTRUE);
     fDiele->GetTrackFilter().AddCuts(trackCuts);
@@ -97,7 +97,7 @@ AliESDtrackCuts *SetupESDtrackCuts()
   esdTrackCuts->SetRequireTPCRefit(kTRUE);
   esdTrackCuts->SetRequireITSRefit(kTRUE);
   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
-  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
+  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
   
   esdTrackCuts->SetMinNClustersTPC(90);
   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
@@ -140,12 +140,12 @@ void InitHistograms()
   histos->UserHistogram("Track","dXY","dXY;dXY [cm];#tracks",500,-1.,1.,AliDielectronVarManager::kImpactParXY);
   histos->UserHistogram("Track","dZ","dZ;dZ [cm];#tracks",600,-3.,3.,AliDielectronVarManager::kImpactParZ);
   histos->UserHistogram("Track","Eta_Phi","Eta Phi Map; Eta; Phi;#tracks",
-                        200,-1,1,200,0,6.285,AliDielectronVarManager::kEta,AliDielectronVarManager::kPhi);
+                        100,-1,1,144,0,6.285,AliDielectronVarManager::kEta,AliDielectronVarManager::kPhi);
   
   histos->UserHistogram("Track","dEdx_P","dEdx;P [GeV];TPC signal (arb units);#tracks",
-                        400,0.2,20.,200,0.,200.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCsignal,kTRUE);
+                        200,0.2,20.,100,0.,200.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCsignal,kTRUE);
   histos->UserHistogram("Track","TPCnSigmaEle_P","TPC number of sigmas Electrons;P [GeV];TPC number of sigmas Electrons;#tracks",
-                        400,0.2,20.,200,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle,kTRUE);
+                        200,0.2,20.,100,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle,kTRUE);
   
   //add histograms to Pair classes
   histos->UserHistogram("Pair","InvMass","Inv.Mass;Inv. Mass [GeV];#pairs",