From 39eb6410a249ed4366e81abe67aa44161f413cc3 Mon Sep 17 00:00:00 2001 From: mfloris Date: Thu, 27 Jan 2011 14:34:23 +0000 Subject: [PATCH] Fitting code from Lee and updated scripts --- PWG2/SPECTRA/LambdaK0PbPb/FitControl.h | 73 ++ PWG2/SPECTRA/LambdaK0PbPb/FitLambdaSpectrum.C | 749 +++++++++--------- PWG2/SPECTRA/LambdaK0PbPb/MultYields2.C | 357 +++++++++ PWG2/SPECTRA/LambdaK0PbPb/PtMassAna2.C | 437 ++++++++++ PWG2/SPECTRA/LambdaK0PbPb/run.sh | 83 +- 5 files changed, 1312 insertions(+), 387 deletions(-) create mode 100644 PWG2/SPECTRA/LambdaK0PbPb/FitControl.h create mode 100644 PWG2/SPECTRA/LambdaK0PbPb/MultYields2.C create mode 100644 PWG2/SPECTRA/LambdaK0PbPb/PtMassAna2.C diff --git a/PWG2/SPECTRA/LambdaK0PbPb/FitControl.h b/PWG2/SPECTRA/LambdaK0PbPb/FitControl.h new file mode 100644 index 00000000000..440eda0b13d --- /dev/null +++ b/PWG2/SPECTRA/LambdaK0PbPb/FitControl.h @@ -0,0 +1,73 @@ +#include "TMath.h" + +class FitControl : public TObject { + protected: + Double_t mPtUpper; // Upper and lower pt limits + Int_t mBinUpper; // Upper and lower bin limits + Int_t mBinLower; + Int_t mPolyOrder; // Polynomial order - 0=constant, 1=linear, 2=quadratic, others not supported + Int_t mRebinFactor; // Rebinning factor + Double_t mMinMass; // Min and max to use as fitting range + Double_t mMaxMass; + + public: + // This data member moved here due to use in sorting functions + Double_t mPtLower; + + FitControl(){ + mPtUpper=0.0; mPtLower=0.0; mPolyOrder = 2; mRebinFactor = 1; + mBinLower=0;mBinUpper=0; + }; + FitControl(Float_t ptLo, Float_t ptUp, Int_t polyO, Int_t RBF){ + mPtUpper=ptUp; mPtLower=ptLo; mPolyOrder = polyO; mRebinFactor = RBF; + mBinLower=0;mBinUpper=0; + }; + FitControl(Double_t ptLo, Double_t ptUp, Int_t polyO, Int_t RBF){ + mPtUpper=ptUp; mPtLower=ptLo; mPolyOrder = polyO; mRebinFactor = RBF; + mMinMass=1.085; mMaxMass=1.17; + mBinLower=0;mBinUpper=0; + }; + FitControl(Double_t ptLo, Double_t ptUp, Int_t polyO, Int_t RBF, Double_t m1, Double_t m2){ + mPtUpper=ptUp; mPtLower=ptLo; mPolyOrder = polyO; mRebinFactor = RBF; + mMinMass=m1; mMaxMass=m2; + mBinLower=0;mBinUpper=0; + }; + FitControl(Int_t BinLo, Int_t BinUp, Int_t polyO, Int_t RBF){ + mPtUpper=0.0; mPtLower=0.0; mPolyOrder = polyO; mRebinFactor = RBF; + mBinLower=BinLo;mBinUpper=BinUp; + }; + ~FitControl(){;}; + + // Functions for sorting + Bool_t IsEqual(const TObject *obj) const {return mPtLower == ((FitControl*)obj)->mPtLower;}; //Not sure whether this one reqd for sorting + Bool_t IsSortable() const { return kTRUE; }; + Int_t Compare(const TObject *obj) const + { + if ( mPtLower < ((FitControl*)obj)->mPtLower) return -1; + if (mPtLower > ((FitControl*)obj)->mPtLower) return 1; + return 0; + }; + + Int_t RebinFactor(){return mRebinFactor;}; + Int_t BinLower(){return mBinLower;}; + Int_t BinUpper(){return mBinUpper;}; + Double_t PtUpper(){return mPtUpper;}; + Double_t PtLower(){return mPtLower;}; + Double_t MinMass(){return mMinMass;}; + Double_t MaxMass(){return mMaxMass;}; + + + Bool_t FixedQuad(){if(mPolyOrder < 2) {return kTRUE;} + else {return kFALSE;}; + }; + Bool_t FixedLin(){if(mPolyOrder < 1) {return kTRUE;} + else {return kFALSE;}; + }; + Double_t DPt(){return mPtUpper-mPtLower;}; + void CalcBinLimits(Int_t BinsPerGeV){ + mBinLower = TMath::Nint(1.+mPtLower*BinsPerGeV); + mBinUpper = TMath::Nint(mPtUpper*BinsPerGeV); + }; + + ClassDef(FitControl,1); +}; diff --git a/PWG2/SPECTRA/LambdaK0PbPb/FitLambdaSpectrum.C b/PWG2/SPECTRA/LambdaK0PbPb/FitLambdaSpectrum.C index 8db83223911..5aca6e3fe1a 100644 --- a/PWG2/SPECTRA/LambdaK0PbPb/FitLambdaSpectrum.C +++ b/PWG2/SPECTRA/LambdaK0PbPb/FitLambdaSpectrum.C @@ -1,123 +1,129 @@ static int myAliceBlue = TColor::GetColor(0,0,192); static int myAliceRed = TColor::GetColor(192,0,0); -void FitLambdaSpectrum(const char* name = "Strange_merge_LHC10h.root",const char* HistName = "h2PtVsMassLambda", const char* SaveName = "SPECTRA_MK_Lambda_7TeV.root") +void FitLambdaSpectrum(const char* name, const char * listName = "clambdak0Histo_00", const char* HistName = "h2PtVsMassLambda", const char* SaveName = "SPECTRA_MK_Lambda_7TeV.root") ////////////////////////zaciatok FitLambdaSpectrum////////////////////////////////////////////////////////// { -gROOT->SetStyle("Plain"); - + gROOT->SetStyle("Plain"); + gROOT->LoadMacro("run.C"); + InitAndLoadLibs(); - TFile *myFile = new TFile(name); - myFile->ls(); - TList *myList = myFile->Get("clist"); + TFile *myFile = new TFile(name); + myFile->ls(); + TList *myList = myFile->Get(listName); // FIXME: find suffix? - TCanvas *myCan = new TCanvas("spectra","lambda"); - myCan->Draw(); - myCan->cd(); + if(!myList) { + cout << "Cannot get list " << listName << " from file " << name << endl; + return; + } - TLatex *myKind = new TLatex(0.15,0.92,"ALICE: LHC10b Pass1 7 TeV"); - myKind->SetNDC(); - myKind->SetTextSize(0.03); - myKind->SetTextColor(myAliceBlue); - myKind->Draw(); + TCanvas *myCan = new TCanvas("spectra","lambda"); + myCan->Draw(); + myCan->cd(); + TLatex *myKind = new TLatex(0.15,0.92,"ALICE: LHC10b Pass1 7 TeV"); // FIXME + myKind->SetNDC(); + myKind->SetTextSize(0.03); + myKind->SetTextColor(myAliceBlue); + myKind->Draw(); - TText *textTitle = new TText(0.6,0.86,""); - textTitle->SetNDC(); - textTitle->SetTextSize(0.04); - textTitle->Draw(); - TPad *myPad1 = new TPad("myPad1","myPad1",0,0,1,0.85); - myPadSetUp(myPad1,0.12,0.02,0.02,0.15); - myPad1->Draw(); - myPad1->cd(); + TText *textTitle = new TText(0.6,0.86,""); + textTitle->SetNDC(); + textTitle->SetTextSize(0.04); + textTitle->Draw(); -// gStyle->SetOptTitle(0); -// gStyle->SetOptStat(0); + TPad *myPad1 = new TPad("myPad1","myPad1",0,0,1,0.85); + myPadSetUp(myPad1,0.12,0.02,0.02,0.15); + myPad1->Draw(); + myPad1->cd(); -TH1F * SaveParameters= new TH1F("SaveParameters","parameters",6,0,6); + // gStyle->SetOptTitle(0); + // gStyle->SetOptStat(0); -///zaciatok binovania///////////// -const Int_t range=33; + TH1F * SaveParameters= new TH1F("SaveParameters","parameters",6,0,6); -Double_t b[range]; + ///zaciatok binovania///////////// + const Int_t range=33; -b[0]=0; -b[1]=0.5; -for (Int_t i = 2;i<17;i++) -{ -b[i]=b[i-1]+0.1; -cout<=1)&&(rBin<16)) {iLoBin = iLoBin+2; iHiBin = iHiBin+2;} + if ((rBin>=1)&&(rBin<16)) {iLoBin = iLoBin+2; iHiBin = iHiBin+2;} -if (rBin==16) {iLoBin = iLoBin+2; iHiBin = iHiBin+4;} + if (rBin==16) {iLoBin = iLoBin+2; iHiBin = iHiBin+4;} -if ((rBin>=17)&&(rBin<26)){iLoBin = iLoBin+4; iHiBin = iHiBin+4;} + if ((rBin>=17)&&(rBin<26)){iLoBin = iLoBin+4; iHiBin = iHiBin+4;} -if (rBin==26) {iLoBin = iLoBin+4; iHiBin = iHiBin+10;} + if (rBin==26) {iLoBin = iLoBin+4; iHiBin = iHiBin+10;} -if (rBin==27) {iLoBin = iLoBin+10; iHiBin = iHiBin+10;} + if (rBin==27) {iLoBin = iLoBin+10; iHiBin = iHiBin+10;} -if (rBin==28) {iLoBin = iLoBin+10; iHiBin = iHiBin+10;} + if (rBin==28) {iLoBin = iLoBin+10; iHiBin = iHiBin+10;} -if (rBin==29) {iLoBin = iLoBin+10; iHiBin = iHiBin+20;} + if (rBin==29) {iLoBin = iLoBin+10; iHiBin = iHiBin+20;} -if (rBin==30) {iLoBin = iLoBin+20; iHiBin = iHiBin+30;} + if (rBin==30) {iLoBin = iLoBin+20; iHiBin = iHiBin+30;} -if (rBin==31) {iLoBin = iLoBin+30; iHiBin = iHiBin+80;} + if (rBin==31) {iLoBin = iLoBin+30; iHiBin = iHiBin+80;} - char saveFileName[60]; - char saveFileName1[60]; - sprintf(saveFileName,"hm_InvariantMassLambdafit_bin%d.gif",rBin); - sprintf(saveFileName1,"hm_InvariantMassLambdafit_bin%d.root",rBin); - myCan->SaveAs(saveFileName); - myCan->SaveAs(saveFileName1); - myCan->Clear(); + char saveFileName[60]; + char saveFileName1[60]; + sprintf(saveFileName,"hm_InvariantMassLambdafit_bin%d.gif",rBin); + sprintf(saveFileName1,"hm_InvariantMassLambdafit_bin%d.root",rBin); + myCan->SaveAs(saveFileName); + myCan->SaveAs(saveFileName1); + myCan->Clear(); - } -//////////////////////////koniec cyklu cez jednotlive biny///////////////////////////// + } + //////////////////////////koniec cyklu cez jednotlive biny///////////////////////////// -//////////////////////////vykreslenie a ulozenie histogramu/////////////////////////// - DrawSpectrumLambda->Draw(); + //////////////////////////vykreslenie a ulozenie histogramu/////////////////////////// + DrawSpectrumLambda->Draw(); - TFile *SPECTRA_MK_Lambda = new TFile(SaveName,"RECREATE"); - SPECTRA_MK_Lambda->WriteObject(DrawSpectrumLambda,"DrawSpectrumLambda"); + TFile *SPECTRA_MK_Lambda = new TFile(SaveName,"RECREATE"); + SPECTRA_MK_Lambda->WriteObject(DrawSpectrumLambda,"DrawSpectrumLambda"); - SPECTRA_MK_Lambda->Close(); -//////////////////////////koniec ulozenia a vykreslenia///////////////////////////////// + SPECTRA_MK_Lambda->Close(); + //////////////////////////koniec ulozenia a vykreslenia///////////////////////////////// } @@ -134,154 +140,154 @@ if (rBin==31) {iLoBin = iLoBin+30; iHiBin = iHiBin+80;} void myBinCounting(TList *myList = 0, TH1F *DrawSpectrumLambda = 0, Int_t rBin = 0, Int_t iLoBin =0, Int_t iHiBin =0, const char* HistName =0,TH1F *SaveParameters = 0){ - TH2F *h2PtVsMassLambda = (TH2F*)myList->FindObject(HistName); - TH1D *h1MassLambdaPtBin = h2PtVsMassLambda->ProjectionX("h1MassLambdaPtBin",iLoBin,iHiBin); - h1MassLambdaPtBin->SetFillColor(0); -// h1MassLambdaPtBin->SetTitle("h1MassLambdaPtBin;invariant mass (GeV/c^{2});Counts"); - //h1MassLambdaPtBin->Draw("SAME"); + TH2F *h2PtVsMassLambda = (TH2F*)myList->FindObject(HistName); + TH1D *h1MassLambdaPtBin = h2PtVsMassLambda->ProjectionX("h1MassLambdaPtBin",iLoBin,iHiBin); + h1MassLambdaPtBin->SetFillColor(0); + // h1MassLambdaPtBin->SetTitle("h1MassLambdaPtBin;invariant mass (GeV/c^{2});Counts"); + //h1MassLambdaPtBin->Draw("SAME"); -///////////////////////////////////////Gauss+2nd pol fit//////////////////////////////////////////// - -TH1D *myHistGauss = h1MassLambdaPtBin; -//Double_t hMax=myHistGauss->GetMaximum(); -//myHistGauss->SetMaximum(hMax); -myHistGauss->Draw(); -// double myLevel = 500, mySlope = 0, mySlope2 = 0 ; -// double myNorm = 100, myMean = 0.497, myWidth = 0.003; - double myRange[2]={0,0}; - - // Gaussain Fit - - double myLevel = 0, mySlope = 0, mySlope2 = 0; + ///////////////////////////////////////Gauss+2nd pol fit//////////////////////////////////////////// + + TH1D *myHistGauss = h1MassLambdaPtBin; + //Double_t hMax=myHistGauss->GetMaximum(); + //myHistGauss->SetMaximum(hMax); + myHistGauss->Draw(); + // double myLevel = 500, mySlope = 0, mySlope2 = 0 ; + // double myNorm = 100, myMean = 0.497, myWidth = 0.003; + double myRange[2]={0,0}; + + // Gaussain Fit + + double myLevel = 0, mySlope = 0, mySlope2 = 0; double myNorm = 0, myMean = 1.116, myWidth = 0.003; -// myNorm=(myHistGauss->GetMaximum()*TMath::Sqrt(2*TMath::Pi())); - myNorm=myHistGauss->GetMaximum(); -// myRange[0]=1.09; myRange[1]=1.14; - myRange[0]=1.095; myRange[1]=1.135; - - if(rBin==4) - { - //myLevel =-2.46E5; mySlope = 3.39E5; mySlope2 = -1.03E5; myNorm = 5.12E3; - myRange[0]=1.10; myRange[1]=1.14; - } - - if(rBin==3) - { - myLevel =-2.91E6; mySlope = 5.21E6; mySlope2 = -2.33E6; myNorm = 9.32E3; - myRange[0]=1.104; myRange[1]=1.14; - } - - char cRange[30],cInterval[30],cIntegration[30],cFitInfo[100], cFitInfo2[100], cFitInfo3[100]; - - TF1 *fPolGauss = new TF1("fpolgauss",myPolGauss,myRange[0],myRange[1],6); - fPolGauss->SetTitle("The gaussian function"); - fPolGauss->SetParNames("level","slope","slope2","norm","mean","width"); - fPolGauss->SetParameters(myLevel,mySlope,mySlope2,myNorm,myMean,myWidth); - fPolGauss->SetLineStyle(2); - fPolGauss->SetLineWidth(2); - fPolGauss->SetLineColor(602); - /// fPolGauss->Draw("SAME"); - myHistGauss->Fit(fPolGauss,"IREM"); - fPolGauss->Draw("SAME"); - - //myHistGauss->Fit(fPolGauss,"LIREM"); - - myNorm = fPolGauss->GetParameter(3); - myMean = fPolGauss->GetParameter(4); - myWidth = fPolGauss->GetParameter(5); - - - -/////store parameters for next fit if own will fall//// -//if ((myMean-(4*myWidth))<(myMean+(4*myWidth))) -//if (myWidth>0.0005) -if ((myWidth>0.0005)&&(myWidth<0.004)) -{ -for (Int_t i=1;i<7;i++) SaveParameters->SetBinContent(i,0); -SaveParameters->SetBinContent(1,fPolGauss->GetParameter(0)); -SaveParameters->SetBinContent(2,fPolGauss->GetParameter(1)); -SaveParameters->SetBinContent(3,fPolGauss->GetParameter(2)); -SaveParameters->SetBinContent(4,fPolGauss->GetParameter(3)); -SaveParameters->SetBinContent(5,fPolGauss->GetParameter(4)); -SaveParameters->SetBinContent(6,fPolGauss->GetParameter(5)); + // myNorm=(myHistGauss->GetMaximum()*TMath::Sqrt(2*TMath::Pi())); + myNorm=myHistGauss->GetMaximum(); + // myRange[0]=1.09; myRange[1]=1.14; + myRange[0]=1.095; myRange[1]=1.135; + + if(rBin==4) + { + //myLevel =-2.46E5; mySlope = 3.39E5; mySlope2 = -1.03E5; myNorm = 5.12E3; + myRange[0]=1.10; myRange[1]=1.14; + } + if(rBin==3) + { + myLevel =-2.91E6; mySlope = 5.21E6; mySlope2 = -2.33E6; myNorm = 9.32E3; + myRange[0]=1.104; myRange[1]=1.14; + } -} -/////////////////////////////////////////////////////// + char cRange[30],cInterval[30],cIntegration[30],cFitInfo[100], cFitInfo2[100], cFitInfo3[100]; + + TF1 *fPolGauss = new TF1("fpolgauss",myPolGauss,myRange[0],myRange[1],6); + fPolGauss->SetTitle("The gaussian function"); + fPolGauss->SetParNames("level","slope","slope2","norm","mean","width"); + fPolGauss->SetParameters(myLevel,mySlope,mySlope2,myNorm,myMean,myWidth); + fPolGauss->SetLineStyle(2); + fPolGauss->SetLineWidth(2); + fPolGauss->SetLineColor(602); + /// fPolGauss->Draw("SAME"); + myHistGauss->Fit(fPolGauss,"IREM"); + fPolGauss->Draw("SAME"); + + //myHistGauss->Fit(fPolGauss,"LIREM"); + + myNorm = fPolGauss->GetParameter(3); + myMean = fPolGauss->GetParameter(4); + myWidth = fPolGauss->GetParameter(5); + + + + /////store parameters for next fit if own will fall//// + //if ((myMean-(4*myWidth))<(myMean+(4*myWidth))) + //if (myWidth>0.0005) + if ((myWidth>0.0005)&&(myWidth<0.004)) + { + for (Int_t i=1;i<7;i++) SaveParameters->SetBinContent(i,0); + SaveParameters->SetBinContent(1,fPolGauss->GetParameter(0)); + SaveParameters->SetBinContent(2,fPolGauss->GetParameter(1)); + SaveParameters->SetBinContent(3,fPolGauss->GetParameter(2)); + SaveParameters->SetBinContent(4,fPolGauss->GetParameter(3)); + SaveParameters->SetBinContent(5,fPolGauss->GetParameter(4)); + SaveParameters->SetBinContent(6,fPolGauss->GetParameter(5)); -///////////if parameters are wrong use previous///// -//if ((myMean-(4*myWidth))>=(myMean+(4*myWidth))) -//if (myWidth<=0.0005) -if ((myWidth<=0.0005)||(myWidth>=0.004)) -{ -cout<<"Previous parameters"<GetBinContent(i)<GetBinContent(1); -mySlope=SaveParameters->GetBinContent(2); -mySlope2=SaveParameters->GetBinContent(3); -myNorm=SaveParameters->GetBinContent(4); -myMean=SaveParameters->GetBinContent(5); -myWidth=SaveParameters->GetBinContent(6); -cout<<"nova sirka"<=(myMean+(4*myWidth))) + //if (myWidth<=0.0005) + if ((myWidth<=0.0005)||(myWidth>=0.004)) + { + cout<<"Previous parameters"<GetBinContent(i)<GetBinContent(1); + mySlope=SaveParameters->GetBinContent(2); + mySlope2=SaveParameters->GetBinContent(3); + myNorm=SaveParameters->GetBinContent(4); + myMean=SaveParameters->GetBinContent(5); + myWidth=SaveParameters->GetBinContent(6); + + cout<<"nova sirka"<SetTitle("The gaussian function"); - fGauss->SetParNames("norm","mean","width"); - fGauss->SetParameters(myNorm,myMean,myWidth); - fGauss->SetLineColor(myAliceRed); - fGauss->SetFillColor(myAliceRed); - fGauss->SetFillStyle(3005); - fGauss->SetLineWidth(1); - fGauss->Draw("SAME"); - - int NumSigmaSignal = 4; - int NumSigmaBackground =6; - double integralWidth = NumSigmaSignal*myWidth; - - double binWidth = myHistGauss->GetBinWidth(1); - double integralGauss = fGauss->Integral(myMean-integralWidth,myMean+integralWidth); - sprintf(cRange,"%d#sigma range (GeV/c^{2})",NumSigmaSignal); - sprintf(cInterval,"[%.3f;%.3f]",myMean-integralWidth,myMean+integralWidth); - sprintf(cIntegration,"integral: %.0f",integralGauss/binWidth); - sprintf(cFitInfo,"#Chi^{2}/ndf = %.1f/%d",fPolGauss->GetChisquare(),fPolGauss->GetNDF()); -// sprintf(cFitInfo2,"mean = %.4f #pm %.4f",fPolGauss->GetParameter(4),fPolGauss->GetParError(4)); -// sprintf(cFitInfo3,"width = %.5f #pm %.5f",fPolGauss->GetParameter(5),fPolGauss->GetParError(5)); + + TF1 *fGauss = new TF1("fgauss",myGauss,myRange[0],myRange[1],3); + fGauss->SetTitle("The gaussian function"); + fGauss->SetParNames("norm","mean","width"); + fGauss->SetParameters(myNorm,myMean,myWidth); + fGauss->SetLineColor(myAliceRed); + fGauss->SetFillColor(myAliceRed); + fGauss->SetFillStyle(3005); + fGauss->SetLineWidth(1); + fGauss->Draw("SAME"); + + int NumSigmaSignal = 4; + int NumSigmaBackground =6; + double integralWidth = NumSigmaSignal*myWidth; + + double binWidth = myHistGauss->GetBinWidth(1); + double integralGauss = fGauss->Integral(myMean-integralWidth,myMean+integralWidth); + sprintf(cRange,"%d#sigma range (GeV/c^{2})",NumSigmaSignal); + sprintf(cInterval,"[%.3f;%.3f]",myMean-integralWidth,myMean+integralWidth); + sprintf(cIntegration,"integral: %.0f",integralGauss/binWidth); + sprintf(cFitInfo,"#Chi^{2}/ndf = %.1f/%d",fPolGauss->GetChisquare(),fPolGauss->GetNDF()); + // sprintf(cFitInfo2,"mean = %.4f #pm %.4f",fPolGauss->GetParameter(4),fPolGauss->GetParError(4)); + // sprintf(cFitInfo3,"width = %.5f #pm %.5f",fPolGauss->GetParameter(5),fPolGauss->GetParError(5)); -cout<GetParameter(4); //stredna hodnota gausovskeho fitu -Double_t GaussSigma =myWidth; //fPolGauss->GetParameter(5); //sigma gausovskeho fitu + Double_t GaussMean =myMean; //fPolGauss->GetParameter(4); //stredna hodnota gausovskeho fitu + Double_t GaussSigma =myWidth; //fPolGauss->GetParameter(5); //sigma gausovskeho fitu -Double_t LeftSide = myHistGauss->FindBin(GaussMean - NumSigmaSignal*GaussSigma); // hranice signalu -Double_t RightSide = myHistGauss->FindBin(GaussMean + NumSigmaSignal*GaussSigma); + Double_t LeftSide = myHistGauss->FindBin(GaussMean - NumSigmaSignal*GaussSigma); // hranice signalu + Double_t RightSide = myHistGauss->FindBin(GaussMean + NumSigmaSignal*GaussSigma); -Double_t NumberBinsSignal = (RightSide - LeftSide) + 1; //pocet binov ktore obsahuju signal + Double_t NumberBinsSignal = (RightSide - LeftSide) + 1; //pocet binov ktore obsahuju signal - Double_t LeftB_Side_Left = myHistGauss->FindBin(myRange[0]); //bacground on left side of peak - Double_t LeftB_Side_Right = myHistGauss->FindBin(GaussMean - NumSigmaBackground*GaussSigma); + Double_t LeftB_Side_Left = myHistGauss->FindBin(myRange[0]); //bacground on left side of peak + Double_t LeftB_Side_Right = myHistGauss->FindBin(GaussMean - NumSigmaBackground*GaussSigma); - Double_t RightB_Side_Left = myHistGauss->FindBin(GaussMean + NumSigmaBackground*GaussSigma); //bacground on Right side of peak - Double_t RightB_Side_Right = myHistGauss->FindBin(myRange[1]); + Double_t RightB_Side_Left = myHistGauss->FindBin(GaussMean + NumSigmaBackground*GaussSigma); //bacground on Right side of peak + Double_t RightB_Side_Right = myHistGauss->FindBin(myRange[1]); -Double_t NumberBinsBackground = (LeftB_Side_Right - LeftB_Side_Left) + 1 + (RightB_Side_Right - RightB_Side_Left) + 1; // pocet binov co + Double_t NumberBinsBackground = (LeftB_Side_Right - LeftB_Side_Left) + 1 + (RightB_Side_Right - RightB_Side_Left) + 1; // pocet binov co // @@ -308,14 +314,14 @@ Double_t NumberBinsBackground = (LeftB_Side_Right - LeftB_Side_Left) + 1 + (Righ nCountNoise +=h1MassLambdaPtBin ->GetBinContent(iBin); myHistNoise->SetBinContent(iBin,h1MassLambdaPtBin->GetBinContent(iBin)); myHistNoise->SetBinError(iBin,h1MassLambdaPtBin->GetBinError(iBin)); - } + } for (Int_t iBin = RightB_Side_Left; iBin <= RightB_Side_Right; iBin++) { nCountNoise +=h1MassLambdaPtBin ->GetBinContent(iBin); myHistNoise->SetBinContent(iBin,h1MassLambdaPtBin->GetBinContent(iBin)); myHistNoise->SetBinError(iBin,h1MassLambdaPtBin->GetBinError(iBin)); - } + } // Create a histo with the Signal in 4 sigma. int totInHistSignal = 0; @@ -329,7 +335,7 @@ Double_t NumberBinsBackground = (LeftB_Side_Right - LeftB_Side_Left) + 1 + (Righ for (Int_t iBin = LeftSide; iBin <= RightSide; iBin++){ myHistSignal->SetBinContent(iBin,h1MassLambdaPtBin->GetBinContent(iBin)); myHistSignal->SetBinError(iBin,h1MassLambdaPtBin->GetBinError(iBin)); - totInHistSignal += h1MassLambdaPtBin->GetBinContent(iBin); + totInHistSignal += h1MassLambdaPtBin->GetBinContent(iBin); } @@ -339,26 +345,26 @@ Double_t NumberBinsBackground = (LeftB_Side_Right - LeftB_Side_Left) + 1 + (Righ // -cout<<"zaciatok background fit casti "<FixParameter(2,(GaussMean - 6*GaussSigma)); + // fitNoise->FixParameter(3,(GaussMean + 6*GaussSigma)); + // } + + + TF1*fitNoise = new TF1("fitNoise",myPol2,myRange[0],myRange[1],5); + fitNoise->SetParLimits(2,-1E8,0); + fitNoise->SetParameter(0,fPolGauss->GetParameter(0)); + fitNoise->SetParameter(1,fPolGauss->GetParameter(1)); + fitNoise->SetParameter(2,fPolGauss->GetParameter(2)); -// if (rOrderPoly == 1 ) { -// TF1*fitNoise = new TF1("fitNoise",myPol1,LeftB_Side_Left,RightB_Side_Right,4); -// -// fitNoise->FixParameter(2,(GaussMean - 6*GaussSigma)); -// fitNoise->FixParameter(3,(GaussMean + 6*GaussSigma)); -// } - - - TF1*fitNoise = new TF1("fitNoise",myPol2,myRange[0],myRange[1],5); - fitNoise->SetParLimits(2,-1E8,0); - fitNoise->SetParameter(0,fPolGauss->GetParameter(0)); - fitNoise->SetParameter(1,fPolGauss->GetParameter(1)); - fitNoise->SetParameter(2,fPolGauss->GetParameter(2)); - - fitNoise->FixParameter(3,(GaussMean - NumSigmaBackground*GaussSigma)); - fitNoise->FixParameter(4,(GaussMean + NumSigmaBackground*GaussSigma)); + fitNoise->FixParameter(3,(GaussMean - NumSigmaBackground*GaussSigma)); + fitNoise->FixParameter(4,(GaussMean + NumSigmaBackground*GaussSigma)); @@ -366,86 +372,86 @@ cout<Fit("fitNoise","irem+"); - //myHistNoise->Fit("fitNoise","Lirem+"); - - Double_t hMax=myHistGauss->GetMaximum(); - myHistNoise->SetMaximum(hMax); - - myHistNoise->GetYaxis()->SetTitle("counts"); - myHistNoise->GetYaxis()->SetTitleOffset(1.2); - myHistNoise->SetMaximum(hMax); - myHistNoise->Draw("HIST"); - myHistSignal->SetFillColor(myAliceRed); - myHistSignal->Draw("A BAR E SAME"); - myHistGauss->SetLineWidth(1); - myHistGauss->Draw("A H E SAME"); - - //fitNoise->Draw("same"); - - //**************************************** - // Bin Counting - //**************************************** - - // Noise Under Peak - -// TF1*fNoise = new TF1("fNoise","[0]+x*[1]",LeftB_Side_Left,RightB_Side_Right); - // fNoise->FixParameter(0,fitNoise->GetParameter(0)); -// fNoise->FixParameter(1,fitNoise->GetParameter(1)); - - - - TF1*fNoise = new TF1("fNoise","[0]+x*[1]+x*x*[2]",myRange[0],myRange[1]); - fNoise->FixParameter(0,fitNoise->GetParameter(0)); - fNoise->FixParameter(1,fitNoise->GetParameter(1)); - fNoise->FixParameter(2,fitNoise->GetParameter(2)); - - - fNoise->Draw("same"); - -// Double_t nNoiseUnderPeak = fNoise->Integral(GaussMean - NumSigmaSignal*GaussSigma,GaussMean + NumSigmaSignal*GaussSigma)/(myHistNoise->GetBinWidth(1)); -Double_t lBinWidth = myHistNoise->GetBinWidth(1); -cout<<"sirka binu"<< lBinWidth<Integral(h1MassLambdaPtBin->GetBinCenter(LeftSide)-0.5*lBinWidth,h1MassLambdaPtBin->GetBinCenter(RightSide)+0.5*lBinWidth)/lBinWidth; - - - Double_t nNoiseUnderPeakErr = TMath::Sqrt(nNoiseUnderPeak); - - - // Signal - Double_t Signal[2]; // Signal[0]=signal and Signal[1]=error - Signal[0] = Signal[1]=0.0; - Signal[0] = totInHistSignal-nNoiseUnderPeak; - Signal[1] = TMath::Sqrt(Signal[0]+nNoiseUnderPeak); - -// DrawSpectrumLambda->SetBinContent(rBin,Signal[0]/DrawSpectrumLambda->GetBinWidth(rBin); -// DrawSpectrumLambda->SetBinError(rBin,Signal[1]/DrawSpectrumLambda->GetBinWidth(rBin)); -DrawSpectrumLambda->SetBinContent(rBin,Signal[0]); -DrawSpectrumLambda->SetBinError(rBin,Signal[1]); - - printf("BoInfo: Noise intervals [%.3f;%.3f],[%.3f;%.3f]\n",myRange[0],(GaussMean - NumSigmaBackground*GaussSigma),(GaussMean +NumSigmaBackground *GaussSigma),myRange[1]); - printf("BoInfo: nBinNoise(left+right)=%d nCountNoise(left+right)=%f \n",NumberBinsBackground,nCountNoise); - printf("BoInfo: nBinCentral=%d totInHistSignal=%f \n",NumberBinsSignal,totInHistSignal); - printf("BoInfo: Noise under peak %.2f +-%.2f \n",nNoiseUnderPeak, nNoiseUnderPeakErr); - printf("BoInfo: Signal %f \n",Signal[0]); - - - // printf only and keep only on plot the Signal, noise and S/N. - Char_t cIntervalSignal[50], cSignalOverNoise[50], cNoiseUnderPeak[50], cFitNoiseInfo[50]; - Double_t SignalOverNoise = 0.0; - Double_t SignalOverNoiseErr = 0.0; - - sprintf(cIntervalSignal,"Signal [%.3f;%.3f] = %.0f #pm %.0f",LeftSide,RightSide,Signal[0], Signal[1]); - sprintf(cNoiseUnderPeak,"Noise under peak = %.0f #pm %.0f",nNoiseUnderPeak, nNoiseUnderPeakErr); - SignalOverNoise = Signal[0]/nNoiseUnderPeak; - SignalOverNoiseErr = TMath::Sqrt(Signal[0]+2*nNoiseUnderPeak); - sprintf(cSignalOverNoise,"S/N= %.2f #pm %.2f",SignalOverNoise,SignalOverNoiseErr); - cout<<"cSignalOverNoise"<GetChisquare(),fitNoise->GetNDF()); + myHistNoise->Fit("fitNoise","irem+"); + //myHistNoise->Fit("fitNoise","Lirem+"); + + Double_t hMax=myHistGauss->GetMaximum(); + myHistNoise->SetMaximum(hMax); + + myHistNoise->GetYaxis()->SetTitle("counts"); + myHistNoise->GetYaxis()->SetTitleOffset(1.2); + myHistNoise->SetMaximum(hMax); + myHistNoise->Draw("HIST"); + myHistSignal->SetFillColor(myAliceRed); + myHistSignal->Draw("A BAR E SAME"); + myHistGauss->SetLineWidth(1); + myHistGauss->Draw("A H E SAME"); + + //fitNoise->Draw("same"); + + //**************************************** + // Bin Counting + //**************************************** + + // Noise Under Peak + + // TF1*fNoise = new TF1("fNoise","[0]+x*[1]",LeftB_Side_Left,RightB_Side_Right); + // fNoise->FixParameter(0,fitNoise->GetParameter(0)); + // fNoise->FixParameter(1,fitNoise->GetParameter(1)); + + + + TF1*fNoise = new TF1("fNoise","[0]+x*[1]+x*x*[2]",myRange[0],myRange[1]); + fNoise->FixParameter(0,fitNoise->GetParameter(0)); + fNoise->FixParameter(1,fitNoise->GetParameter(1)); + fNoise->FixParameter(2,fitNoise->GetParameter(2)); + + + fNoise->Draw("same"); + + // Double_t nNoiseUnderPeak = fNoise->Integral(GaussMean - NumSigmaSignal*GaussSigma,GaussMean + NumSigmaSignal*GaussSigma)/(myHistNoise->GetBinWidth(1)); + Double_t lBinWidth = myHistNoise->GetBinWidth(1); + cout<<"sirka binu"<< lBinWidth<Integral(h1MassLambdaPtBin->GetBinCenter(LeftSide)-0.5*lBinWidth,h1MassLambdaPtBin->GetBinCenter(RightSide)+0.5*lBinWidth)/lBinWidth; + + + Double_t nNoiseUnderPeakErr = TMath::Sqrt(nNoiseUnderPeak); + + + // Signal + Double_t Signal[2]; // Signal[0]=signal and Signal[1]=error + Signal[0] = Signal[1]=0.0; + Signal[0] = totInHistSignal-nNoiseUnderPeak; + Signal[1] = TMath::Sqrt(Signal[0]+nNoiseUnderPeak); + + // DrawSpectrumLambda->SetBinContent(rBin,Signal[0]/DrawSpectrumLambda->GetBinWidth(rBin); + // DrawSpectrumLambda->SetBinError(rBin,Signal[1]/DrawSpectrumLambda->GetBinWidth(rBin)); + DrawSpectrumLambda->SetBinContent(rBin,Signal[0]); + DrawSpectrumLambda->SetBinError(rBin,Signal[1]); + + printf("BoInfo: Noise intervals [%.3f;%.3f],[%.3f;%.3f]\n",myRange[0],(GaussMean - NumSigmaBackground*GaussSigma),(GaussMean +NumSigmaBackground *GaussSigma),myRange[1]); + printf("BoInfo: nBinNoise(left+right)=%d nCountNoise(left+right)=%f \n",NumberBinsBackground,nCountNoise); + printf("BoInfo: nBinCentral=%d totInHistSignal=%f \n",NumberBinsSignal,totInHistSignal); + printf("BoInfo: Noise under peak %.2f +-%.2f \n",nNoiseUnderPeak, nNoiseUnderPeakErr); + printf("BoInfo: Signal %f \n",Signal[0]); + -cout<GetChisquare(),fitNoise->GetNDF()); + + cout<= abscMin && absc < abscMax) { - TF1::RejectPoint(); - return 0; - } - - double ordo = level + slope * absc; - return ordo; - -} - -double myPol2(double *x, double *pars) { - - double absc = x[0]; - double level = pars[0]; - double slope = pars[1]; - double slope2 = pars[2]; - double abscMin = pars[3]; - double abscMax = pars[4]; - - if (absc >= abscMin && absc < abscMax) { - TF1::RejectPoint(); - return 0; - } - - double ordo = level + slope * absc +slope2 * absc*absc; - return ordo; - -} +double myGauss(double *x, double *par){ + double absc = x[0]; + double norm = par[0]; + double mean = par[1]; + double width = par[2]; + // double mygauss = (norm/TMath::Sqrt(2*TMath::Pi()))*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width)); + double mygauss = (norm)*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width)); + return mygauss; +} + +double myPolGauss(double *x, double *pars){ + double absc = x[0]; + double level = pars[0]; + double slope = pars[1]; + double slope2 = pars[2]; + double norm = pars[3]; + double mean = pars[4]; + double width = pars[5]; + // double ordo = level + slope * absc +slope2 * absc*absc + (norm/TMath::Sqrt(2*TMath::Pi()))*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width)); + double ordo = level + slope * absc +slope2 * absc*absc + (norm)*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width)); + return ordo; +} + +double myPol1(double *x, double *pars) { + + double absc = x[0]; + double level = pars[0]; + double slope = pars[1]; + double abscMin = pars[2]; + double abscMax = pars[3]; + + if (absc >= abscMin && absc < abscMax) { + TF1::RejectPoint(); + return 0; + } + + double ordo = level + slope * absc; + return ordo; + +} + +double myPol2(double *x, double *pars) { + + double absc = x[0]; + double level = pars[0]; + double slope = pars[1]; + double slope2 = pars[2]; + double abscMin = pars[3]; + double abscMax = pars[4]; + + if (absc >= abscMin && absc < abscMax) { + TF1::RejectPoint(); + return 0; + } + + double ordo = level + slope * absc +slope2 * absc*absc; + return ordo; + +} +void LoadLibs(){ + gSystem->Load("libCore.so"); + gSystem->Load("libTree.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libPhysics.so"); + gSystem->Load("libSTEERBase"); + gSystem->Load("libESD"); + gSystem->Load("libAOD"); + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); +} diff --git a/PWG2/SPECTRA/LambdaK0PbPb/MultYields2.C b/PWG2/SPECTRA/LambdaK0PbPb/MultYields2.C new file mode 100644 index 00000000000..b688753cd1c --- /dev/null +++ b/PWG2/SPECTRA/LambdaK0PbPb/MultYields2.C @@ -0,0 +1,357 @@ +#include "FitControl.h" + +void MultYields2(TH3F *PtMassMult, Int_t particleMode, Int_t MultBin, Char_t* label){ + + TString *tLabel = new TString(label); + // Do .L on command line first + //gROOT->LoadMacro("macros/PtMassAna2.C"); +/* Modifications to produce a single uncorrected spectrum in a particular mult. bin for K0 or Lambda +Old ratio code preserved in MultYieldsRatio.C +Dec 2003 */ + Char_t* part; //for name of particle + Float_t BR; //branching ratio - check them + if (particleMode==0){ + part = "K0"; + BR=0.686; + } else if (particleMode==1){ + part = "Lambda"; + BR=0.639; + } else if (particleMode==2){ + part = "AntiLambda"; + BR=0.639; + } else if (particleMode==3){ + part = "Lambda+antiLambda"; + BR=0.639; + } else if (particleMode==4){ + part = "Xi"; + BR=1.0; // Should be Lam ->p pi Br + } else if (particleMode ==6){ + part = "Omega"; + BR=1.0; // Should be Lam->p pi * Om -> K Lam + } + + +// //Setup. Multiplicity bin boundaries, N_binary scaling numbers +// //These are latest numbers from M. Miller Aug03 +// Float_t NBin[4] = {7.5,3.99,10.2,15.04}; +// Float_t Veff[4] = {0.925,0.875,0.999,1.0}; +// // These are dAuFilt numbers +// // Float_t Nev1=5.03505e+06, Nev2=2.4501e+06, Nev3=2.23543e+06, NevMB=9.722229e+06; +// Float_t MultLo[4]={0,0,10,17}; +// Float_t MultHi[4]={120,9,16,120}; +// //These are dAuFilt2 numbers (discrepancy with above, changed GoodEvent BUT that shouldn't decrease Nevent) +// //New N events figures from the 25cm cut file +// if (tLabel->Contains("25cm")) { +// Float_t Nev[4]={4.586559e+06,2.37952e+06,1.11382e+06,1.09322e+06}; +// } +// else { +// Float_t Nev[4]={9.590239e+06,5.0683e+06,2.31685e+06,2.20508e+06}; +// } + +// Redo the possible multiplicity bin stuff for Cu+Cu SVT test +// As a start just use one bin + Float_t NBin[1] = {63.5}; // Not really that useful for this study + // used 0-80% bin from Johan's page + Float_t Veff[1] = {0.91}; // TBD from data? + Float_t Nev[1] = {34712502.0}; // For MB +// Float_t Nev[1] = {2727356.0}; // For HM + Float_t MultLo[1] = {0}; + Float_t MultHi[1] = {300}; + +TString title[1]={"MinimumBias"}; + //Make 2D projections from the 3D histogram + //Minbias (i.e. everything) + PtMassMult->GetZaxis()->SetRange(MultLo[MultBin],MultHi[MultBin]); + TH2F* hPtMass = (TH2F*)PtMassMult->Project3D("XY");// FIX:MF + hPtMass->SetTitle("PtMass"); + + hPtMass->Draw(); + + TObjArray *controllerArray = new TObjArray(40,1); // 2nd arg means can count from 1! + + //Here probably need switch-case, depending on mult bin and particle + // LoPt, HiPt, polynomial order, rebinning factor + + ///// LAMBDA and LAMBDA+ANTILAMBDA (Combination) + + if(particleMode == 1 || particleMode == 3){ // Lambda or Lambda+Anti-Lambda + // controllerArray->AddLast(new FitControl(0.1,0.2, 1,2)); + // controllerArray->AddLast(new FitControl(0.2,0.3, 1,2)); + // controllerArray->AddLast(new FitControl(0.3,0.4, 1,1, 1.096,1.16)) + // controllerArray->AddLast(new FitControl(0.4,0.5, 2,1, 1.095,1.15)); + if(tLabel->Contains("SVT")){ + controllerArray->AddLast(new FitControl(0.4,0.5, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(0.5,0.6, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(0.6,0.7, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(0.7,0.8, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(0.8,1.0, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(1.0,1.2, 2,1)); + controllerArray->AddLast(new FitControl(1.2,1.4, 2,1)); + controllerArray->AddLast(new FitControl(1.4,1.6, 2,1)); + controllerArray->AddLast(new FitControl(1.6,1.8, 2,1)); + controllerArray->AddLast(new FitControl(1.8,2.0, 2,2)); + controllerArray->AddLast(new FitControl(2.0,2.2, 2,2, 1.095,1.17)); + controllerArray->AddLast(new FitControl(2.2,2.4, 2,2, 1.095,1.16)); +// controllerArray->AddLast(new FitControl(2.4,2.6, 2,1, 1.095,1.17)); +// controllerArray->AddLast(new FitControl(2.6,2.8, 1,2, 1.095,1.17)); +// controllerArray->AddLast(new FitControl(2.8,3.0, 1,2, 1.095,1.17)); +// controllerArray->AddLast(new FitControl(3.0,3.5, 1,2, 1.098,1.16)); +// controllerArray->AddLast(new FitControl(3.5,4.0, 1,2, 1.095,1.17)); + }else{ + controllerArray->AddLast(new FitControl(0.5,0.6, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(0.6,0.7, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(0.7,0.8, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(0.8,1.0, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(1.0,1.2, 2,1)); + controllerArray->AddLast(new FitControl(1.2,1.4, 2,1)); + controllerArray->AddLast(new FitControl(1.4,1.6, 2,1)); + controllerArray->AddLast(new FitControl(1.6,1.8, 2,1)); + controllerArray->AddLast(new FitControl(1.8,2.0, 2,1)); + controllerArray->AddLast(new FitControl(2.0,2.2, 2,1, 1.095,1.17)); + controllerArray->AddLast(new FitControl(2.2,2.4, 2,1, 1.095,1.16)); + controllerArray->AddLast(new FitControl(2.4,2.6, 2,1, 1.095,1.17)); + controllerArray->AddLast(new FitControl(2.6,2.8, 1,2, 1.095,1.17)); + controllerArray->AddLast(new FitControl(2.8,3.0, 1,2, 1.095,1.17)); + controllerArray->AddLast(new FitControl(3.0,3.5, 1,2, 1.098,1.16)); + controllerArray->AddLast(new FitControl(3.5,4.0, 1,2, 1.095,1.17)); + } +// if (MultBin==2){controllerArray->AddLast(new FitControl(4.0,4.5, 1,2, 1.090,1.17));} +// else if (MultBin==3){controllerArray->AddLast(new FitControl(4.0,4.5, 0,3, 1.090,1.17));} +// else {controllerArray->AddLast(new FitControl(4.0,4.5, 1,2, 1.095,1.17));} +// controllerArray->AddLast(new FitControl(4.5,5.0, 0,4, 1.095,1.17)); +// if (MultBin==1){ controllerArray->AddLast(new FitControl(5.0,5.5, 0,5, 1.095,1.17));} +// else if (MultBin==2){controllerArray->AddLast(new FitControl(5.0,5.5, 0,6, 1.085,1.16));} +// else {controllerArray->AddLast(new FitControl(5.0,5.5, 0,4, 1.095,1.17));} +// if(MultBin==3){ +// controllerArray->AddLast(new FitControl(5.5,6.0, 0,5, 1.085,1.165)); +// } else if(MultBin==2){ +// controllerArray->AddLast(new FitControl(5.5,6.0, 0,5, 1.085,1.17)); +// } else{ +// controllerArray->AddLast(new FitControl(5.5,6.0, 0,6, 1.085,1.17)); +// } + // controllerArray->AddLast(new FitControl(6.0,6.5, 0,6)); + // controllerArray->AddLast(new FitControl(6.5,7.0, 0,10)); + // The above, bins from 500 MeV to 6 GeV works for CENTRAL Lambdas + // i.e all chisq/dof for backgrounds ~<2, no wild fits, perhaps sig/bkgrd too + // high for 2.8-3.0 GeV bin + } + /// ANTI LAMBDA ----> + else if (particleMode == 2){ // Anti-Lambdas + controllerArray->AddLast(new FitControl(0.3,0.4, 2,1, 1.095,1.15)); + if(MultBin==3) { controllerArray->AddLast(new FitControl(0.4,0.5, 2,1, 1.095,1.16));} + else { controllerArray->AddLast(new FitControl(0.4,0.5, 2,1, 1.095,1.15));} + controllerArray->AddLast(new FitControl(0.5,0.6, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(0.6,0.7, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(0.7,0.8, 2,1, 1.090,1.15)); + controllerArray->AddLast(new FitControl(0.8,1.0, 2,1, 1.095,1.15)); + controllerArray->AddLast(new FitControl(1.0,1.2, 2,1, 1.090,1.17)); + + if (MultBin==1){ controllerArray->AddLast(new FitControl(1.2,1.4, 2,1, 1.090,1.16));} + else{controllerArray->AddLast(new FitControl(1.2,1.4, 2,1, 1.090,1.17));} + controllerArray->AddLast(new FitControl(1.4,1.6, 2,1, 1.090,1.17)); + controllerArray->AddLast(new FitControl(1.6,1.8, 2,1, 1.090,1.17)); + controllerArray->AddLast(new FitControl(1.8,2.0, 2,1, 1.090,1.17)); + controllerArray->AddLast(new FitControl(2.0,2.2, 2,1, 1.095,1.17)); + controllerArray->AddLast(new FitControl(2.2,2.4, 2,1, 1.095,1.16)); + if(MultBin==1){ controllerArray->AddLast(new FitControl(2.4,2.6, 1,1, 1.095,1.17));} + else if(MultBin==3){ controllerArray->AddLast(new FitControl(2.4,2.6, 1,1, 1.095,1.17));} + else { controllerArray->AddLast(new FitControl(2.4,2.6, 2,1, 1.095,1.17));} + controllerArray->AddLast(new FitControl(2.6,2.8, 1,2, 1.095,1.17)); + if(MultBin==2){ controllerArray->AddLast(new FitControl(2.8,3.0, 1,2, 1.090,1.17));} + else if(MultBin==3){ controllerArray->AddLast(new FitControl(2.8,3.0, 1,2, 1.090,1.17));} + else{ controllerArray->AddLast(new FitControl(2.8,3.0, 1,2, 1.095,1.17));} + controllerArray->AddLast(new FitControl(3.0,3.5, 1,2, 1.095,1.17)); + controllerArray->AddLast(new FitControl(3.5,4.0, 1,2, 1.095,1.17)); + controllerArray->AddLast(new FitControl(4.0,4.5, 1,2, 1.090,1.17)); + if (MultBin==1){ + // controllerArray->AddLast(new FitControl(4.5,5.0, 0,5, 1.095,1.17)); + controllerArray->AddLast(new FitControl(4.5,5.0, 0,3, 1.085,1.17)); + controllerArray->AddLast(new FitControl(5.0,5.5, 0,3, 1.090,1.17)); + //controllerArray->AddLast(new FitControl(5.5,6.0, 0,6, 1.085,1.17)); + } else if (MultBin==2){ + controllerArray->AddLast(new FitControl(4.5,5.0, 0,4, 1.095,1.17)); + controllerArray->AddLast(new FitControl(5.0,5.5, 0,4, 1.095,1.17)); + controllerArray->AddLast(new FitControl(5.5,6.0, 0,4, 1.085,1.17)); + } else if (MultBin==3){ + controllerArray->AddLast(new FitControl(4.5,5.0, 0,4, 1.095,1.17)); + controllerArray->AddLast(new FitControl(5.0,5.5, 0,4, 1.095,1.17)); + controllerArray->AddLast(new FitControl(5.5,6.0, 0,4, 1.085,1.18)); + } else{ + controllerArray->AddLast(new FitControl(4.5,5.0, 0,4, 1.095,1.17)); + controllerArray->AddLast(new FitControl(5.0,5.5, 0,4, 1.095,1.17)); + controllerArray->AddLast(new FitControl(5.5,6.0, 0,6, 1.085,1.17)); + } + } // end if anti-Lambda + else if (particleMode == 0){ // K0s case + if(tLabel->Contains("SVT")){ + controllerArray->AddLast(new FitControl(0.4,0.5, 1,1, 0.441,0.558)); + controllerArray->AddLast(new FitControl(0.5,0.6, 1,1, 0.442,0.559)); + controllerArray->AddLast(new FitControl(0.6,0.7, 1,1, 0.442,0.558)); + controllerArray->AddLast(new FitControl(0.7,0.8, 1,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(0.8,1.0, 1,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(1.0,1.2, 2,1, 0.442,0.558)); + controllerArray->AddLast(new FitControl(1.2,1.4, 2,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(1.4,1.6, 2,1, 0.442,0.558)); + controllerArray->AddLast(new FitControl(1.6,1.8, 2,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(1.8,2.0, 1,2, 0.44,0.56)); + controllerArray->AddLast(new FitControl(2.0,2.2, 1,2, 0.44,0.56)); + controllerArray->AddLast(new FitControl(2.2,2.4, 1,2, 0.44,0.56)); + controllerArray->AddLast(new FitControl(2.4,2.6, 1,2, 0.44,0.56)); + controllerArray->AddLast(new FitControl(2.6,2.8, 1,4, 0.44,0.56)); + controllerArray->AddLast(new FitControl(2.8,3.0, 1,4, 0.44,0.56)); + }else{ + //controllerArray->AddLast(new FitControl(0.1,0.2, 1,2, 0.44,0.56)); + //controllerArray->AddLast(new FitControl(0.2,0.3, 1,2, 0.44,0.56)); + //controllerArray->AddLast(new FitControl(0.3,0.4, 1,1, 0.442,0.559)); + controllerArray->AddLast(new FitControl(0.4,0.5, 1,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(0.5,0.6, 1,1, 0.442,0.559)); + controllerArray->AddLast(new FitControl(0.6,0.7, 1,1, 0.442,0.558)); + controllerArray->AddLast(new FitControl(0.7,0.8, 1,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(0.8,1.0, 1,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(1.0,1.2, 2,1, 0.442,0.558)); + controllerArray->AddLast(new FitControl(1.2,1.4, 2,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(1.4,1.6, 2,1, 0.442,0.558)); + controllerArray->AddLast(new FitControl(1.6,1.8, 2,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(1.8,2.0, 1,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(2.0,2.2, 1,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(2.2,2.4, 1,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(2.4,2.6, 1,1, 0.44,0.56)); + controllerArray->AddLast(new FitControl(2.6,2.8, 1,2, 0.44,0.56)); + controllerArray->AddLast(new FitControl(2.8,3.0, 1,2, 0.44,0.56)); + controllerArray->AddLast(new FitControl(3.0,3.5, 1,2, 0.44,0.56)); + controllerArray->AddLast(new FitControl(3.5,4.0, 1,2, 0.44,0.56)); + // controllerArray->AddLast(new FitControl(4.0,4.5, 1,2, 0.418,0.578)); + // controllerArray->AddLast(new FitControl(4.5,5.0, 0,2, 0.418,0.578)); + // controllerArray->AddLast(new FitControl(5.0,5.5, 0,4, 0.418,0.578)); + // controllerArray->AddLast(new FitControl(5.5,6.0, 0,4, 0.418,0.578)); + // controllerArray->AddLast(new FitControl(5.5,6.0, 0,4, 0.418,0.568)); + // controllerArray->AddLast(new FitControl(6.0,6.5, 0,4, 0.418,0.578)); + + } + // Above values all good for PERIPHERAL K0 + // For MINBIAS K0 there is a problem for the 1.4-1.6 GeV bin - this bin has possibly the most asymmetric background + // For MIDCENTRAL K0, bins 1.4-1.6 and 1.6-1.8 GeV have problems also bin 4.5-5 GeV should probably be rebinned (peak looks split) + // For CENTRAL K0s bin 1.4-1.6 fails. The two bins before that have S/N ratios which deviate from the trend. + //Reduced fit range for 1.4-1.6 bin means working OK for CENTRAL now. + + // controllerArray->AddLast(new FitControl(6.5,7.0, 0,10)); + } else if (particleMode == 4) { //Xi case + //controllerArray->AddLast(new FitControl(0.5,0.7, 1,1, 1.28,1.45)); //signal not visible with + controllerArray->AddLast(new FitControl(0.7,0.9, 1,1, 1.285,1.45)); + controllerArray->AddLast(new FitControl(0.9,1.1, 1,1, 1.285,1.45)); + controllerArray->AddLast(new FitControl(1.1,1.3, 1,1, 1.27,1.45)); + controllerArray->AddLast(new FitControl(1.3,1.5, 1,1, 1.27,1.45)); + controllerArray->AddLast(new FitControl(1.5,1.7, 1,1, 1.27,1.45)); + controllerArray->AddLast(new FitControl(1.7,2.0, 1,1, 1.27,1.45)); + controllerArray->AddLast(new FitControl(2.0,2.5, 1,1, 1.27,1.44)); + controllerArray->AddLast(new FitControl(2.5,3.0, 1,1, 1.27,1.45)); + controllerArray->AddLast(new FitControl(3.0,3.5, 1,1, 1.27,1.44)); + controllerArray->AddLast(new FitControl(3.5,4.0, 1,1, 1.27,1.45)); + controllerArray->AddLast(new FitControl(4.0,4.5, 2,1, 1.27,1.44)); + controllerArray->AddLast(new FitControl(4.5,5.0, 2,1, 1.27,1.45)); + } + TString *tLabel = new TString(label); + // Special setup for doing the comparison to dEdx kaons (different binning) + if(tLabel->Contains("KdEdxComp")){ + controllerArray->Delete(); // removes all that went before! + controllerArray->AddLast(new FitControl(0.25,0.35,1,2,0.418,0.578)); + controllerArray->AddLast(new FitControl(0.35,0.45,1,2,0.418,0.578)); + controllerArray->AddLast(new FitControl(0.45,0.55,1,2,0.418,0.578)); + controllerArray->AddLast(new FitControl(0.55,0.65,1,2,0.418,0.568)); + controllerArray->AddLast(new FitControl(0.65,0.75,2,2,0.418,0.578)); + controllerArray->AddLast(new FitControl(0.75,0.85,2,2,0.418,0.578)); + controllerArray->AddLast(new FitControl(0.85,0.95,2,2,0.418,0.568)); + } + // Special setup for doing comparison to Hai's K0s (different binning) + if(tLabel->Contains("HJBins")){ + controllerArray->Delete(); // removes all that went before! + controllerArray->AddLast(new FitControl(0.4,0.6,1,1,0.418,0.578)); + controllerArray->AddLast(new FitControl(0.6,0.8,1,1,0.418,0.568)); + controllerArray->AddLast(new FitControl(0.8,1.0,1,1,0.418,0.578)); + controllerArray->AddLast(new FitControl(1.0,1.2,2,1,0.418,0.578)); + controllerArray->AddLast(new FitControl(1.2,1.4,2,1,0.418,0.578)); + controllerArray->AddLast(new FitControl(1.4,1.6,2,1,0.418,0.568)); + controllerArray->AddLast(new FitControl(1.6,1.8,2,1,0.418,0.578)); + controllerArray->AddLast(new FitControl(1.8,2.0,2,1,0.418,0.568)); + controllerArray->AddLast(new FitControl(2.0,2.5,2,2,0.418,0.558)); + controllerArray->AddLast(new FitControl(2.5,3.0,1,2,0.418,0.568)); + controllerArray->AddLast(new FitControl(3.0,3.5,1,2,0.418,0.578)); + controllerArray->AddLast(new FitControl(3.5,4.5,1,2,0.418,0.578)); + controllerArray->AddLast(new FitControl(4.5,6.0,0,4,0.418,0.578)); + } + controllerArray->Sort(); + + // Make the proper label to pass in for use in labelling the diagnostics + // saved canvas files and histos + Char_t fulllabel[80]; + sprintf(fulllabel,"%s%s",title[MultBin].Data(),label); + + //Slice up projection into mass histograms to extract yield +TH1F* hYield = (TH1F*)PtMassAna2(hPtMass,particleMode,controllerArray->GetEntries(),controllerArray,fulllabel); + +// CORRECTIONS - only do number of events for now +hYield->Scale(1.0/Nev[MultBin]); //Divide by the number of events +//hYield->Scale(Veff[MultBin]); //Multiply by the vertex efficiency effectively increaing number of events + //(since Veff<1) therefore decreases yield +//hYield->Scale(1.0/BR); //Divide by branching ratio (again increases yield since BR<1) + //hYield->Scale(1.0/(2*TMath::Pi())); // Always plot 1/2pi ... + +Char_t yieldTitle[80]; +sprintf(yieldTitle,"Uncorrected %s yield: %s",part,title[MultBin].Data()); +hYield->SetTitle(yieldTitle); +hYield->SetXTitle("p_{t} / [GeV/c]"); +hYield->SetYTitle("1/Nev.dN/dp_{t}"); + +// Create plots + +Char_t fileNameBase[80]; +sprintf(fileNameBase,"Masses%s_%s%s",part,title[MultBin].Data(),label); +Char_t fileNamePng[80]; +sprintf(fileNamePng,"%s.png",fileNameBase); +Char_t fileNameEps[80]; +sprintf(fileNameEps,"%s.eps",fileNameBase); +Char_t fileNamePdf[80]; +sprintf(fileNamePdf,"%s.pdf",fileNameBase); + +c1->SaveAs(fileNamePng); +c1->SaveAs(fileNameEps); +c1->SaveAs(fileNamePdf); + +//c1->Clear(); + +//c1->SetLogy(); +TCanvas *cYield = new TCanvas("Yield","Corrected Yield",600,400); +cYield->cd(); +//cYield->SetLogy(); +hYield->SetStats(kFALSE); +hYield->Draw(); + cYield->SetLogy(); + cYield->Update(); +// hRC_MB->SetMarkerStyle(20); +// hRC_MB->SetMarkerColor(4); +// hRC_MB->Scale(NBinMB/NBin3); +Char_t fnametext[80]; +sprintf(fnametext,"Yield%s_%s%s",part,title[MultBin].Data(),label); +Char_t fnamePng[80]; +sprintf(fnamePng,"%s.png",fnametext); +c1->SaveAs(fnamePng); +Char_t fnameEps[80]; +sprintf(fnameEps,"%s.eps",fnametext); +c1->SaveAs(fnameEps); + +TH1F* hScYield = hYield->Clone("ScYield"); +Char_t scYieldTitle[80]; +sprintf(scYieldTitle," Scaled %s",hYield->GetTitle()); +hScYield->SetTitle(scYieldTitle); +//SCALING for scaled yield only divide by mean Nbin (scaled yield is therefore smaller) +hScYield->Scale(1/NBin[MultBin]); + +Char_t fnameRoot[80]; +sprintf(fnameRoot,"%s.root",fnametext); +TFile *YieldFile = new TFile(fnameRoot,"RECREATE"); +hYield->Write(); +hScYield->Write(); +YieldFile->Close(); + + controllerArray->Delete(); +} diff --git a/PWG2/SPECTRA/LambdaK0PbPb/PtMassAna2.C b/PWG2/SPECTRA/LambdaK0PbPb/PtMassAna2.C new file mode 100644 index 00000000000..77e105124ec --- /dev/null +++ b/PWG2/SPECTRA/LambdaK0PbPb/PtMassAna2.C @@ -0,0 +1,437 @@ +#include "FitControl.h" + +// Function for fitting back ground away from peak +Float_t quad(Double_t *x, Double_t *par){ + if (x[0] > par[3] && x[0] < par[4]){ + TF1::RejectPoint(); + return 0; + } + return par[0] + par[1]*x[0] + par[2]*x[0]*x[0]; //a+bx+cx**2 +} + +TH1F* PtMassAna2(TH2F *PtMass, Int_t mode, const Int_t NControl, TObjArray *ControlArray,Char_t* fulllabel){ + + //TString *tLabel = new TString(fulllabel); //not reqd + + //Arguments TH2F *PtMass - 2D histogram of some variable versus mass + // Int_t mode - mode for switching things depending on particle + // const Int_t NControl - number of bins to project into which must + // be the same as the number of controller objects (see next arg). + // TObjArray *ContolArray - holds the controllers objects, 1 per bin + + //const Int_t NDefault=27; // Used in default array (change as required) + // Trap user errors + if (ControlArray->GetEntries() != NControl){ + cout << "PtMassAna2: Problem: Specified number of bins and number of bins in supplied TObjArray do not match" << endl; + return 0; + } + //Set up mode mapping to particle and names for histograms + if(mode == 0){ + const Char_t* part = "K0"; + const Char_t* partName = "K^{0}_{s}"; + } + else if (mode == 1 ){ + const Char_t* part = "LAM"; + const Char_t* partName = "#Lambda"; + } + else if (mode == 2 ){//since lambda and anti-Lambda have same limits etc. + const Char_t* part = "LAM"; + const Char_t* partName = "#bar{#Lambda}"; + } + else if (mode ==3) {// This is for combined Lambda + anti-Lambda + const Char_t* part = "LAM"; + const Char_t* partName = "#Lambda + #bar{#Lambda}"; + } + else if (mode ==4) { + const Char_t* part = "XI"; + const Char_t* partName = "#Xi"; + } + else{ + cout << "Mode not recognized, defaulting to K0" << endl; + const Char_t* part = "K0"; + } + cout << "Debug info. Particle: " << part << endl; + + Float_t xxlo, xxhi; //Exclusion limits for fitting + // FUTURE - these can also be part of FitControl to allow different regions + // in different bins + Float_t NSigmaEx=4.0; // Number of sigma to exclude from bkgd fit + Float_t NSigmaInt=3.5; // Number of sigma to integrate histo over + + // Fitting function for initial combined peak + background fit + TF1 *gausQuad = new TF1("gausQuad","(([1]/([2]*pow(6.2831,0.5)))*[6]*exp(-0.5*pow((x-[0])/[2],2)))+ [3] + [4]*x + [5]*x*x",0.3,1.8); + //>>>>>>>>> NB Need to set proper range when it is used <<<<<<<<<<<<<<<<<<<<< + //>>>>>>>>> AND use 'R' in fit option <<<<<<<<<<<<<<<<<<<<< + gausQuad->SetLineWidth(1); + gausQuad->SetLineColor(kRed); + + // Fitting function for background + TF1* quad = new TF1("quad",quad,0.3,1.8,5); + quad->SetLineWidth(2); + quad->SetLineColor(kBlue); + quad->SetLineStyle(2); + + // Function for background under peak - to be integrated and drawn + TF1* qback = new TF1("qback","[0]+[1]*x+[2]*x*x",0.3,1.8); + qback->SetLineWidth(2); + qback->SetLineColor(kRed); + qback->SetLineStyle(2); + qback->SetFillStyle(3013); + qback->SetFillColor(kRed); + + // Set ranges and limits according to particle + Float_t xmin,xmax; + Float_t defMean,defArea,defWidth,defa0,defa1,defa2,defBW; //defaults + Float_t defMnLo, defMnHi; + if (part=="K0"){ + // FUTURE: xmin, xmax define the range considered in fit so could go into + // FitControl + //xmin=0.418; + //xmax=0.578; + // xxlo=0.46;xxhi=0.53; + // gausQuad->SetParameters(0.4977,2000,0.003,1,1,0,0.001*RBFact); + //gausQuad->SetParameters(0.4977,2000,0.003,1,1,0,0.001); + defMean=0.4977; defArea=1000; defWidth=0.01; + defa0=2; defa1=2; defa2=0; + //gausQuad->SetParLimits(0,0.483,0.523); // Constrain mean + defMnLo = 0.483; defMnHi=0.523; + } + if (part=="LAM"){ + // FUTURE: See above + //xmin=1.085; + //xmax=1.17; + // xxlo=1.10;xxhi=1.135; + // gausQuad->SetParameters(1.115,2000,0.0028,10,100,-100,0.001*RBFact); + // gausQuad->SetParameters(1.115,2000,0.0028,10,100,-100,0.001); + defMean=1.115; defArea=1000; defWidth=0.08; + defa0=100; defa1=-10; defa2=0; + // gausQuad->SetParLimits(0,1.113,1.123); // Constrain mean + defMnLo=1.113; defMnHi=1.123; + } + if (part=="XI"){ + defMean=1.32171; defArea=500; defWidth=0.08; + defa0=100; defa1=-10; defa2=0; + defMnLo=1.301; defMnHi=1.341; + } + //gausQuad->SetRange(xmin,xmax); + gausQuad->SetParLimits(1,0.0,1E7); // Constrain area 0 to 10 million + gausQuad->SetParLimits(2,0.0001,0.04); // Constrain width 0.1 MeV to 40 MeV + gausQuad->SetParNames("Mean","Area1","Sigma1","p0","p1","p2","BinWid"); + + // Deciding how to divide up the canvas + c1->Clear(); + c1->SetWindowSize(1024,768); + Int_t NDraw = NControl; + if(NDraw==16){c1->Divide(4,4);} + elseif(NDraw<=25 && NDraw>20){ c1->Divide(5,5);} + elseif(NDraw<=20 && NDraw>16){ c1->Divide(5,4);} + elseif(NDraw<=15 && NDraw>12){ c1->Divide(5,3);} + elseif(NDraw<=12 && NDraw>8){ c1->Divide(4,3);} + elseif(NDraw<=8 && NDraw>6){ c1->Divide(4,2);} + elseif(NDraw<=6 && NDraw>4){c1->Divide(3,2);} + elseif(NDraw<=4 && NDraw>2){c1->Divide(2,2);} + elseif(NDraw==2){c1->Divide(2,1);} + elseif(NDraw==1){/*do nothing*/;} + else{c1->Divide(6,5);} + + // + // Project out the histograms from 2D into 1D 'slices' + // + Char_t id[20], title[80]; + //cout << "NControl: " << NControl << endl; + TH1D* hMassSlice[NControl]; + + //Arrays to store various quantities which can later be histogrammed. + // const Int_t NBinsArrays = NBins+1-NFirst; // array of 1D projection histograms is counted from 0 but other arrays are used for histograms contents and therefore N+1 are need because element zero is left empty + const Int_t NBinsArrays = NControl+1; + Stat_t sigBkgd[NBinsArrays], errSigBkgd[NBinsArrays]; // Signal+background and error on this + Stat_t bkgd[NBinsArrays], errBkgd[NBinsArrays]; // Background and error on this + Stat_t signal[NBinsArrays], errSignal[NBinsArrays]; // Net signal and error + + Stat_t chi2PerDOF1[NBinsArrays], errChi2PerDOF1[NBinsArrays]; // Chi-squared per defree of freedom from the initial fit + Stat_t chi2PerDOF2[NBinsArrays], errChi2PerDOF2[NBinsArrays]; // Chi-squared per defree of freedom from the 2nd (background) fit + + // REDO won't need this? + // Zero the signal - avoids problems with plotting + for(Int_t j=0;jCalcBinLimits(10); + //Had to introduce this Nint fn otherwise got inconsistencies after type implicit conversion + // BinLo=TMath::Nint(1.+BinPtEdges[N]*20.); //cout << "BinLo: " << BinLo << ", " << 1+BinPtEdges[N]*20. << endl; + // BinHi=TMath::Nint(BinPtEdges[N+1]*20.); //cout << "BinHi: " << BinHi << ", " << BinPtEdges[N+1]*20. << endl; + BinLo=controller->BinLower(); + BinHi=controller->BinUpper(); + Int_t N=ControlArray->IndexOf(controller) - 1; + sprintf(id,"Mass%d",N); + cout << "Mass histo:" << N << " Firstbin: " << BinLo << " Last bin:" << BinHi << " " << controller->PtLower() << "-" << controller->PtUpper() << " GeV" << endl; + //cout << "About to create mass projection " << N << endl; + hMassSlice[N] = PtMass->ProjectionY(id,BinLo,BinHi); + //cout << "Mass projection " << N << " created." << endl; + sprintf(title,"%s Mass, %.2f < p_{t} < %.2f",partName,controller->PtLower(),controller->PtUpper()); + hMassSlice[N]->SetTitle(title); + hMassSlice[N]->SetXTitle("GeV/c^{2}"); + // cout << "Mass projection " << N << " titles set." << endl; + hMassSlice[N]->Rebin(controller->RebinFactor()); + // } // End of creating the projections + Int_t massBinLo,massBinHi; //For mass histogram bins + Int_t NArray; + // Do all the fitting + //REDO this is all one loop (iteration over TObjArray) now + //for(Int_t N=0;Ncd(N+1); //Canvas subdivisions numbered from 1 and NOT 0 + + // Set range - can be different each time + xmin=controller->MinMass(); + xmax=controller->MaxMass(); + gausQuad->SetRange(xmin,xmax); + // Everytime need to set some sensible parameters for both fits + gausQuad->SetParameters(defMean,defArea,defWidth,defa0,defa1,defa2,0.1); + quad->SetParameters(defa0,defa1,defa2,0.49,0.51);// Last 2 parameters are + // the range to exclude but these are fixed later by FixParameter(3,.. etc. + gausQuad->FixParameter(6,hMassSlice[N]->GetBinWidth(1)); + // Release Linear parameter + gausQuad->ReleaseParameter(4); + quad->ReleaseParameter(1); + // Release Quadratic parameter + gausQuad->ReleaseParameter(5); + quad->ReleaseParameter(2); + gausQuad->SetParLimits(0,defMnLo,defMnHi); + gausQuad->SetParLimits(1,0.0,1E7); // Constrain area 0 to 10 million + gausQuad->SetParLimits(2,0.0001,0.04); // Constrain width 0.1 MeV to 40 MeV +// gausQuad->SetParLimits(4,0,0); // Releases linear parameter +// quad->SetParLimits(1,0,0); // Releases linear parameter +// gausQuad->SetParLimits(5,0,0); // Releases quadratic parameter +// quad->SetParLimits(2,0,0); // Releases quadratic parameter + if(controller->FixedLin()){ + quad->FixParameter(1,0.0); + gausQuad->FixParameter(4,0.0); + cout << "Info: linear part fixed to zero." << endl; + gausQuad->SetParLimits(3,0.0,1000); // Ensure constant is +ve + // NB documentation indicates that limits only used when option B + // is specified. Also how are they can they be freed later? + } + if(controller->FixedQuad()){ + quad->FixParameter(2,0.0); + gausQuad->FixParameter(5,0.0); + cout << "Info: quadratic part fixed to zero." << endl; + } + hMassSlice[N]->Fit(gausQuad,"LR"); + hMassSlice[N]->Draw("E,SAME"); + gausQuad->DrawCopy("SAME"); + + //Fit the background: + xxlo=gausQuad->GetParameter(0)-NSigmaEx*gausQuad->GetParameter(2); + xxhi=gausQuad->GetParameter(0)+NSigmaEx*gausQuad->GetParameter(2); + // Limits have to be adjusted to match bins + massBinLo=hMassSlice[N]->FindBin(xxlo); + xxlo = hMassSlice[N]->GetBinLowEdge(massBinLo); + massBinHi=hMassSlice[N]->FindBin(xxhi)+1; + xxhi = hMassSlice[N]->GetBinLowEdge(massBinHi); + quad->SetParameters(gausQuad->GetParameter(3),gausQuad->GetParameter(4),gausQuad->GetParameter(5),xxlo,xxhi); + quad->FixParameter(3,xxlo); + quad->FixParameter(4,xxhi); + quad->SetRange(xmin,xmax); + hMassSlice[N]->Fit(quad,"LRN+"); + quad->SetFillColor(30); + quad->SetFillStyle(3013); + quad->SetRange(xmin,xxlo); + quad->DrawCopy("SAME"); + quad->SetRange(xxhi,xmax); + quad->DrawCopy("SAME"); + + // Draw the background - use a new function since other one not defined everywhere + xxlo=gausQuad->GetParameter(0)-NSigmaInt*gausQuad->GetParameter(2); + xxhi=gausQuad->GetParameter(0)+NSigmaInt*gausQuad->GetParameter(2); + massBinLo=hMassSlice[N]->FindBin(xxlo); + xxlo = hMassSlice[N]->GetBinLowEdge(massBinLo); + massBinHi=hMassSlice[N]->FindBin(xxhi)+1; + xxhi = hMassSlice[N]->GetBinLowEdge(massBinHi); + qback->SetRange(xxlo,xxhi); + qback->SetParameter(0,quad->GetParameter(0)); + qback->SetParameter(1,quad->GetParameter(1)); + qback->SetParameter(2,quad->GetParameter(2)); + qback->DrawCopy("SAME"); + + //Integrate the signal+background (i.e. original histo) + sigBkgd[NArray]=hMassSlice[N]->Integral(massBinLo,massBinHi); + errSigBkgd[NArray]=TMath::Sqrt(sigBkgd[NArray]); + //Integrate background - better to do this analytically ? + bkgd[NArray]=qback->Integral(xxlo,xxhi)/gausQuad->GetParameter(6);//Divide by bin width + errBkgd[NArray]=TMath::Sqrt(bkgd[NArray]); //Get errors from fit parameters? + + // Fill arrays for diagnostic histograms + mean[NArray]=gausQuad->GetParameter(0); + errMean[NArray]=gausQuad->GetParError(0); + area[NArray]=gausQuad->GetParameter(1); + errArea[NArray]=gausQuad->GetParError(1); + sigma[NArray]=gausQuad->GetParameter(2); + errSigma[NArray]=gausQuad->GetParError(2); + chi2PerDOF1[NArray]=gausQuad->GetChisquare()/(Double_t)gausQuad->GetNDF(); + errChi2PerDOF1[NArray]= 0.0; // Don't know how to calc. error for this? + chi2PerDOF2[NArray]=quad->GetChisquare()/(Double_t)quad->GetNDF(); + errChi2PerDOF2[NArray]= 0.0; // Don't know how to calc. error for this? + + BinPtEdges[N]=controller->PtLower(); + // BinPtEdges(N+1)=controller->PtUpper(); + } + BinPtEdges[NBinsArrays-1]=((FitControl*)ControlArray->Last())->PtUpper(); +// for (Int_t jj=0;jjReset("ICE"); + hYields->SetXTitle("p_{t} / (GeV/c)"); + //hYields->SetContent(sigBkgd); + //hYields->SetError(errSigBkgd); + //c1->cd(4); + //hYields->Draw(); + + // Fill the diagnostic histograms: clone, set title, contents, errors. + TH1F* hMeans = hYields->Clone("Means"); + hMeans->SetTitle("Mean from Gaussian Fit"); + hMeans->SetContent(mean); + hMeans->SetError(errMean); + //hMeans->Sumw2(); + TH1F* hSigmas = hYields->Clone("Sigmas"); + hSigmas->SetTitle("Sigma from Gaussian Fit"); + hSigmas->SetContent(sigma); + hSigmas->SetError(errSigma); + TH1F* hAreas = hYields->Clone("Areas"); + hAreas->SetTitle("Area from Gaussian Fit"); + hAreas->SetContent(area); + hAreas->SetError(errArea); + TH1F* hSigBkgd = hYields->Clone("SigBkgd"); + hSigBkgd->SetTitle("Signal + Background (Histo. Integral)"); + hSigBkgd->SetContent(sigBkgd); + hSigBkgd->SetError(errSigBkgd); + TH1F* hBkgd = hYields->Clone("Bkgd"); + hBkgd->SetTitle("Background (Fn. integral)"); + hBkgd->SetContent(bkgd); + hBkgd->SetError(errBkgd); + + hYields->Sumw2(); + hYields->Add(hSigBkgd,hBkgd,1.0,-1.0); + + TH1F* hSoverB = hYields->Clone("SoverB"); + hSoverB->SetTitle("Signal to Background ratio"); + hSoverB->Sumw2(); + hSoverB->Divide(hBkgd); + + TH1F* hSoverA = hYields->Clone("SoverA"); // Signal over area + hSoverA->SetTitle("Ratio of Signal to Area from Fit"); + hSoverA->Sumw2(); + hSoverA->Divide(hAreas); + + TH1F* hChi2PerDOF1 = hYields->Clone("Chi2PerDOF1"); + hChi2PerDOF1->SetTitle("Chi-squared per d.o.f. - initial fit"); + hChi2PerDOF1->SetContent(chi2PerDOF1); + hChi2PerDOF1->SetError(errChi2PerDOF1); + + TH1F* hChi2PerDOF2 = hYields->Clone("Chi2PerDOF2"); + hChi2PerDOF2->SetTitle("Chi-squared per d.o.f. - background fit"); + hChi2PerDOF2->SetContent(chi2PerDOF2); + hChi2PerDOF2->SetError(errChi2PerDOF2); + + // Draw the diagnostic histograms on their own canvas + TCanvas* cDiag = new TCanvas("Diag","Diagnostics",600,600); + cDiag->Divide(2,3); + + cDiag->cd(1); + Float_t bookmass; + if(part=="LAM"){ + hMeans->SetMinimum(1.11); + hMeans->SetMaximum(1.13); + bookmass=1.11563; + } else if (part=="K0"){ + hMeans->SetMaximum(0.51); + hMeans->SetMinimum(0.475); + bookmass=0.4976; + } else if (part=="XI") { + hMeans->SetMaximum(1.34); + hMeans->SetMinimum(1.30); + bookmass=1.32171; + } + TLine PDGmass; + PDGmass.SetLineStyle(2); + hMeans->Draw(); + PDGmass.DrawLine(0,bookmass,5,bookmass); + + cDiag->cd(2); + if (part=="K0"){ + hSigmas->SetMaximum(0.012); + }else{ + hSigmas->SetMaximum(0.006);} + hSigmas->SetMinimum(0.0); + hSigmas->Draw(); + + cDiag->cd(3); + hSoverB->SetMinimum(0.0); + hSoverB->Draw(); + + cDiag->cd(4); + hSoverA->SetMinimum(0.5); + hSoverA->SetMaximum(1.5); + hSoverA->Draw(); + + cDiag->cd(5); + hChi2PerDOF1->SetMinimum(0); + hChi2PerDOF1->SetMaximum(6); + hChi2PerDOF1->Draw(); + + cDiag->cd(6); + hChi2PerDOF2->SetMinimum(0); + hChi2PerDOF2->SetMaximum(6); + hChi2PerDOF2->Draw(); + + Char_t fileNameBase[80]; + sprintf(fileNameBase,"Diagnostics%s_%s",part,fulllabel); + Char_t fileNamePng[80]; + sprintf(fileNamePng,"%s.png",fileNameBase); + Char_t fileNameEps[80]; + sprintf(fileNameEps,"%s.eps",fileNameBase); + Char_t fileNamePdf[80]; + sprintf(fileNamePdf,"%s.pdf",fileNameBase); + + cDiag->SaveAs(fileNamePng); + cDiag->SaveAs(fileNameEps); + cDiag->SaveAs(fileNamePdf); + + // Scale result by the bin size to get dN/dpT and by bin centre to get 1/pt... + Float_t y, ey; + for(Int_t j=0;j<=hYields->GetNbinsX();j++){ + y = hYields->GetBinContent(j)/hYields->GetBinWidth(j); + ey = hYields->GetBinError(j)/hYields->GetBinWidth(j); + hYields->SetBinContent(j,y); + hYields->SetBinError(j,ey); + y = hYields->GetBinContent(j)/hYields->GetBinCenter(j); + ey = hYields->GetBinError(j)/hYields->GetBinCenter(j); + hYields->SetBinContent(j,y); + hYields->SetBinError(j,ey); + } + return hYields; + +} diff --git a/PWG2/SPECTRA/LambdaK0PbPb/run.sh b/PWG2/SPECTRA/LambdaK0PbPb/run.sh index dfedcd0a992..c2fd5f7022e 100755 --- a/PWG2/SPECTRA/LambdaK0PbPb/run.sh +++ b/PWG2/SPECTRA/LambdaK0PbPb/run.sh @@ -12,6 +12,12 @@ offset=0 debug=kTRUE option="SAVE" #FIXME:set option suffix="" +fitFolder="LHC10h_000137161_p1_5plus" +fitBin="00" +partID=1 +task=no +fit=no + give_help() { @@ -19,7 +25,7 @@ cat < Run the task Modes (compulsory): 0 local @@ -34,6 +40,7 @@ Available options: (one entry per line). If you use this option you don't need the -d. in the case you are running on CAF, they must have the same path + -x Add extra suffix to files Grid only options -g Plugin Mode [default=$mode] -p Reconstruction pass [default=$pass] @@ -41,14 +48,34 @@ Available options: -p Data set path -n Number of events -w Number of workers [default=$workers] + + * Fit the results * + -f Run the fitting macro in the subfolder of ./output + -b Centrality bin index [default=$fitBin] + -p Fit particle defined by particleID [default=$partID] + 0=K0 + 1=Lambda + 2=Anti-Lambda + 3=Lambda + Anti-Lambda, + 4=Csi + 6=Omega + -x Add extra suffix to files ENDOFGUIDE } -while getopts "r:hd:mg:p:n:w:t:l:" opt; do +while getopts "r:hd:mg:p:n:w:t:l:f:b:x:" opt; do case $opt in r) runMode=$OPTARG + task=yes + ;; + f) + fitFolder=$OPTARG + fit=yes + ;; + b) + fitBin=`printf %2.2d $OPTARG` ;; d) run=$OPTARG @@ -73,6 +100,10 @@ while getopts "r:hd:mg:p:n:w:t:l:" opt; do ;; p) pass=$OPTARG + partID=$OPTARG + ;; + x) + suffix=$OPTARG ;; h) give_help @@ -92,28 +123,38 @@ while getopts "r:hd:mg:p:n:w:t:l:" opt; do done -runlist=$run -if [ "$listfile" != "" ] +if [ "$task" = "yes" ] then - runlist="" - while read line - do - runlist="$runlist $line" - done < $listfile - -fi - -echo "Run list: $runlist" + runlist=$run + if [ "$listfile" != "" ] + then + runlist="" + while read line + do + runlist="$runlist $line" + done < $listfile + fi + echo "Run list: $runlist" + -if [ "$runMode" = "2" ] -then - echo root $ROPT run.C\(\"$run\",\"$pass\",$nev,$offset,$debug,$runMode,$mc,$option,$suffix,$workers,\"$mode\"\) - root $ROPT run.C\(\"$run\",\"$pass\",$nev,$offset,$debug,$runMode,$mc,\"$option\",\"$suffix\",$workers,\"$mode\"\) -else - for run in $runlist - do + if [ "$runMode" = "2" ] + then echo root $ROPT run.C\(\"$run\",\"$pass\",$nev,$offset,$debug,$runMode,$mc,$option,$suffix,$workers,\"$mode\"\) root $ROPT run.C\(\"$run\",\"$pass\",$nev,$offset,$debug,$runMode,$mc,\"$option\",\"$suffix\",$workers,\"$mode\"\) - done + else + for run in $runlist + do + echo root $ROPT run.C\(\"$run\",\"$pass\",$nev,$offset,$debug,$runMode,$mc,$option,$suffix,$workers,\"$mode\"\) + root $ROPT run.C\(\"$run\",\"$pass\",$nev,$offset,$debug,$runMode,$mc,\"$option\",\"$suffix\",$workers,\"$mode\"\) + done + fi +elif [ "$fit" = "yes" ] +then + root FitSpectrum.C\(\"./output/$fitFolder/lambdak0_${fitBin}.root\",\"clambdak0Histo_${fitBin}\",\"$suffix\",$partID\) +else + give_help fi + + + -- 2.39.3