--- /dev/null
+#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);
+};
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<<b[i]<<endl;
-}
+ Double_t b[range];
-for (Int_t i = 17;i<27;i++)
-{
-b[i]=b[i-1]+0.2;
-cout<<b[i]<<endl;
-}
-b[27]=4.5;
-b[28]=5.0;
-b[29]=5.5;
-b[30]=6.5;
-b[31]=8.0;
-b[32]=12.0;
+ b[0]=0;
+ b[1]=0.5;
+ for (Int_t i = 2;i<17;i++)
+ {
+ b[i]=b[i-1]+0.1;
+ cout<<b[i]<<endl;
+ }
+
+ for (Int_t i = 17;i<27;i++)
+ {
+ b[i]=b[i-1]+0.2;
+ cout<<b[i]<<endl;
+ }
+ b[27]=4.5;
+ b[28]=5.0;
+ b[29]=5.5;
+ b[30]=6.5;
+ b[31]=8.0;
+ b[32]=12.0;
-//koniec binovania
-TH1F *DrawSpectrumLambda = new TH1F("DrawSpectrumLambda","#Lambda Spectrum;p_{t} [GeV/c]; N",32,b);
+ //koniec binovania
+ TH1F *DrawSpectrumLambda = new TH1F("DrawSpectrumLambda","#Lambda Spectrum;p_{t} [GeV/c]; N",32,b);
-///////////////////////////////////////////////run over all bins//////////////////////////////////////////////////////////////
- int iLoBin = 15, iHiBin = 16;
-// int iLoBin = 71, iHiBin = 72, hMax = 250;
+ ///////////////////////////////////////////////run over all bins//////////////////////////////////////////////////////////////
+ int iLoBin = 15, iHiBin = 16;
+ // int iLoBin = 71, iHiBin = 72, hMax = 250;
- for (Int_t rBin = 4;rBin<33;rBin++){
-//for (Int_t rBin = 36;rBin<51;rBin++){
+ for (Int_t rBin = 4;rBin<33;rBin++){
+ //for (Int_t rBin = 36;rBin<51;rBin++){
- myBinCounting(myList, DrawSpectrumLambda, rBin,iLoBin ,iHiBin,HistName,SaveParameters);
+ myBinCounting(myList, DrawSpectrumLambda, rBin,iLoBin ,iHiBin,HistName,SaveParameters);
-if ((rBin>=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/////////////////////////////////
}
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 ;\r
-// double myNorm = 100, myMean = 0.497, myWidth = 0.003;\r
- double myRange[2]={0,0}; \r
- \r
- // Gaussain Fit\r
-\r
- double myLevel = 0, mySlope = 0, mySlope2 = 0; \r
+ ///////////////////////////////////////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;\r
- 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;
- }\r
-\r
- char cRange[30],cInterval[30],cIntegration[30],cFitInfo[100], cFitInfo2[100], cFitInfo3[100];\r
- \r
- TF1 *fPolGauss = new TF1("fpolgauss",myPolGauss,myRange[0],myRange[1],6);\r
- fPolGauss->SetTitle("The gaussian function");\r
- fPolGauss->SetParNames("level","slope","slope2","norm","mean","width");\r
- fPolGauss->SetParameters(myLevel,mySlope,mySlope2,myNorm,myMean,myWidth);\r
- fPolGauss->SetLineStyle(2);\r
- fPolGauss->SetLineWidth(2);\r
- fPolGauss->SetLineColor(602);\r
- /// fPolGauss->Draw("SAME");\r
- myHistGauss->Fit(fPolGauss,"IREM");
- fPolGauss->Draw("SAME");
-\r
- //myHistGauss->Fit(fPolGauss,"LIREM");\r
- \r
- myNorm = fPolGauss->GetParameter(3);\r
- myMean = fPolGauss->GetParameter(4);\r
- 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"<<endl;
-for (Int_t i=1;i<7;i++) cout<<SaveParameters->GetBinContent(i)<<endl;
-myLevel=SaveParameters->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"<<myWidth<<endl;
-}\r
-////////////////////////////////////////////////////
+ }
+ ///////////////////////////////////////////////////////
+
+ ///////////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"<<endl;
+ for (Int_t i=1;i<7;i++) cout<<SaveParameters->GetBinContent(i)<<endl;
+ myLevel=SaveParameters->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"<<myWidth<<endl;
+ }
+ ////////////////////////////////////////////////////
- \r
- TF1 *fGauss = new TF1("fgauss",myGauss,myRange[0],myRange[1],3);\r
- fGauss->SetTitle("The gaussian function");\r
- fGauss->SetParNames("norm","mean","width");\r
- fGauss->SetParameters(myNorm,myMean,myWidth);\r
- fGauss->SetLineColor(myAliceRed);\r
- fGauss->SetFillColor(myAliceRed);\r
- fGauss->SetFillStyle(3005);\r
- fGauss->SetLineWidth(1);\r
- fGauss->Draw("SAME");\r
- \r
- int NumSigmaSignal = 4;
- int NumSigmaBackground =6;\r
- double integralWidth = NumSigmaSignal*myWidth;\r
- \r
- double binWidth = myHistGauss->GetBinWidth(1);\r
- double integralGauss = fGauss->Integral(myMean-integralWidth,myMean+integralWidth);\r
- sprintf(cRange,"%d#sigma range (GeV/c^{2})",NumSigmaSignal);\r
- sprintf(cInterval,"[%.3f;%.3f]",myMean-integralWidth,myMean+integralWidth);\r
- sprintf(cIntegration,"integral: %.0f",integralGauss/binWidth);\r
- sprintf(cFitInfo,"#Chi^{2}/ndf = %.1f/%d",fPolGauss->GetChisquare(),fPolGauss->GetNDF());\r
-// sprintf(cFitInfo2,"mean = %.4f #pm %.4f",fPolGauss->GetParameter(4),fPolGauss->GetParError(4));\r
-// sprintf(cFitInfo3,"width = %.5f #pm %.5f",fPolGauss->GetParameter(5),fPolGauss->GetParError(5));\r
+
+ 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<<cRange<<endl;
-cout<<cInterval<<endl;
-cout<<cIntegration<<endl;
-cout<<cFitInfo<<endl;
-//cout<<cFitInfo2<<endl;
-//cout<<cFitInfo3<<endl;
+ cout<<cRange<<endl;
+ cout<<cInterval<<endl;
+ cout<<cIntegration<<endl;
+ cout<<cFitInfo<<endl;
+ //cout<<cFitInfo2<<endl;
+ //cout<<cFitInfo3<<endl;
-//break;
+ //break;
-////////////////////////////////////////koniec Gauss+2nd pol fit////////////////////////////////////
+ ////////////////////////////////////////koniec Gauss+2nd pol fit////////////////////////////////////
-Double_t GaussMean =myMean; //fPolGauss->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
//
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;
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);
}
//
-cout<<"zaciatok background fit casti "<<endl;
-cout<<endl;
-cout<<endl;
+ cout<<"zaciatok background fit casti "<<endl;
+ cout<<endl;
+ cout<<endl;
+
+ // 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));
-// if (rOrderPoly == 1 ) {\r
-// TF1*fitNoise = new TF1("fitNoise",myPol1,LeftB_Side_Left,RightB_Side_Right,4);\r
-//\r
-// fitNoise->FixParameter(2,(GaussMean - 6*GaussSigma));\r
-// fitNoise->FixParameter(3,(GaussMean + 6*GaussSigma)); \r
-// }\r
-\r
-\r
- TF1*fitNoise = new TF1("fitNoise",myPol2,myRange[0],myRange[1],5);\r
- fitNoise->SetParLimits(2,-1E8,0);\r
- fitNoise->SetParameter(0,fPolGauss->GetParameter(0));\r
- fitNoise->SetParameter(1,fPolGauss->GetParameter(1));\r
- fitNoise->SetParameter(2,fPolGauss->GetParameter(2));\r
-\r
- fitNoise->FixParameter(3,(GaussMean - NumSigmaBackground*GaussSigma));\r
- fitNoise->FixParameter(4,(GaussMean + NumSigmaBackground*GaussSigma)); \r
+ fitNoise->FixParameter(3,(GaussMean - NumSigmaBackground*GaussSigma));
+ fitNoise->FixParameter(4,(GaussMean + NumSigmaBackground*GaussSigma));
Double_t inteGrapmyHistSignal = 0;
- myHistNoise->Fit("fitNoise","irem+");\r
- //myHistNoise->Fit("fitNoise","Lirem+");
-\r
- Double_t hMax=myHistGauss->GetMaximum();
- myHistNoise->SetMaximum(hMax);\r
-
- myHistNoise->GetYaxis()->SetTitle("counts");\r
- myHistNoise->GetYaxis()->SetTitleOffset(1.2);
- myHistNoise->SetMaximum(hMax);\r
- myHistNoise->Draw("HIST");\r
- myHistSignal->SetFillColor(myAliceRed);\r
- myHistSignal->Draw("A BAR E SAME");\r
- myHistGauss->SetLineWidth(1);\r
- myHistGauss->Draw("A H E SAME");
-\r
- //fitNoise->Draw("same");\r
-\r
- //****************************************\r
- // Bin Counting\r
- //****************************************\r
-\r
- // Noise Under Peak \r
-\r
-// TF1*fNoise = new TF1("fNoise","[0]+x*[1]",LeftB_Side_Left,RightB_Side_Right); \r
- // fNoise->FixParameter(0,fitNoise->GetParameter(0));\r
-// fNoise->FixParameter(1,fitNoise->GetParameter(1));\r
-\r
-\r
-\r
- TF1*fNoise = new TF1("fNoise","[0]+x*[1]+x*x*[2]",myRange[0],myRange[1]); \r
- fNoise->FixParameter(0,fitNoise->GetParameter(0));\r
- fNoise->FixParameter(1,fitNoise->GetParameter(1));\r
- fNoise->FixParameter(2,fitNoise->GetParameter(2));\r
-\r
-\r
- fNoise->Draw("same");\r
-\r
-// Double_t nNoiseUnderPeak = fNoise->Integral(GaussMean - NumSigmaSignal*GaussSigma,GaussMean + NumSigmaSignal*GaussSigma)/(myHistNoise->GetBinWidth(1));
-Double_t lBinWidth = myHistNoise->GetBinWidth(1);
-cout<<"sirka binu"<< lBinWidth<<endl;
-Double_t nNoiseUnderPeak = fNoise->Integral(h1MassLambdaPtBin->GetBinCenter(LeftSide)-0.5*lBinWidth,h1MassLambdaPtBin->GetBinCenter(RightSide)+0.5*lBinWidth)/lBinWidth;
-
-\r
- Double_t nNoiseUnderPeakErr = TMath::Sqrt(nNoiseUnderPeak);\r
-\r
-
- // Signal\r
- Double_t Signal[2]; // Signal[0]=signal and Signal[1]=error\r
- Signal[0] = Signal[1]=0.0;\r
- Signal[0] = totInHistSignal-nNoiseUnderPeak;\r
- Signal[1] = TMath::Sqrt(Signal[0]+nNoiseUnderPeak);\r
-
-// 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]);
- \r
- printf("BoInfo: Noise intervals [%.3f;%.3f],[%.3f;%.3f]\n",myRange[0],(GaussMean - NumSigmaBackground*GaussSigma),(GaussMean +NumSigmaBackground *GaussSigma),myRange[1]);\r
- printf("BoInfo: nBinNoise(left+right)=%d nCountNoise(left+right)=%f \n",NumberBinsBackground,nCountNoise);\r
- printf("BoInfo: nBinCentral=%d totInHistSignal=%f \n",NumberBinsSignal,totInHistSignal);\r
- printf("BoInfo: Noise under peak %.2f +-%.2f \n",nNoiseUnderPeak, nNoiseUnderPeakErr);\r
- printf("BoInfo: Signal %f \n",Signal[0]);\r
- \r
- \r
- // printf only and keep only on plot the Signal, noise and S/N.\r
- Char_t cIntervalSignal[50], cSignalOverNoise[50], cNoiseUnderPeak[50], cFitNoiseInfo[50];\r
- Double_t SignalOverNoise = 0.0;\r
- Double_t SignalOverNoiseErr = 0.0;\r
- \r
- sprintf(cIntervalSignal,"Signal [%.3f;%.3f] = %.0f #pm %.0f",LeftSide,RightSide,Signal[0], Signal[1]);\r
- sprintf(cNoiseUnderPeak,"Noise under peak = %.0f #pm %.0f",nNoiseUnderPeak, nNoiseUnderPeakErr);\r
- SignalOverNoise = Signal[0]/nNoiseUnderPeak;\r
- SignalOverNoiseErr = TMath::Sqrt(Signal[0]+2*nNoiseUnderPeak);\r
- sprintf(cSignalOverNoise,"S/N= %.2f #pm %.2f",SignalOverNoise,SignalOverNoiseErr);\r
- cout<<"cSignalOverNoise"<<cSignalOverNoise<<endl;\r
- sprintf(cFitNoiseInfo,"#Chi^{2}/ndf = %.1f/%d",fitNoise->GetChisquare(),fitNoise->GetNDF());\r
+ 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<<endl;
+ Double_t nNoiseUnderPeak = fNoise->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<<cIntervalSignal<<endl;
-cout<<cNoiseUnderPeak<<endl;
-cout<<cFitNoiseInfo<<endl;\r
+ // 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"<<cSignalOverNoise<<endl;
+ sprintf(cFitNoiseInfo,"#Chi^{2}/ndf = %.1f/%d",fitNoise->GetChisquare(),fitNoise->GetNDF());
+
+ cout<<cIntervalSignal<<endl;
+ cout<<cNoiseUnderPeak<<endl;
+ cout<<cFitNoiseInfo<<endl;
}
return;
}
-\r
-double myGauss(double *x, double *par){\r
- double absc = x[0];\r
- double norm = par[0];\r
- double mean = par[1];\r
- double width = par[2];\r
-// double mygauss = (norm/TMath::Sqrt(2*TMath::Pi()))*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width));\r
-double mygauss = (norm)*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width));
- return mygauss;\r
-}\r
-\r
-double myPolGauss(double *x, double *pars){\r
- double absc = x[0];\r
- double level = pars[0];\r
- double slope = pars[1];\r
- double slope2 = pars[2];\r
- double norm = pars[3];\r
- double mean = pars[4];\r
- double width = pars[5];\r
-// double ordo = level + slope * absc +slope2 * absc*absc + (norm/TMath::Sqrt(2*TMath::Pi()))*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width));\r
-double ordo = level + slope * absc +slope2 * absc*absc + (norm)*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width));
- return ordo;\r
-}\r
-\r
-double myPol1(double *x, double *pars) {\r
-\r
- double absc = x[0];\r
- double level = pars[0];\r
- double slope = pars[1];\r
- double abscMin = pars[2];\r
- double abscMax = pars[3];\r
-\r
- if (absc >= abscMin && absc < abscMax) {\r
- TF1::RejectPoint();\r
- return 0;\r
- }\r
- \r
- double ordo = level + slope * absc;\r
- return ordo;\r
-\r
-}\r
-\r
-double myPol2(double *x, double *pars) {\r
-\r
- double absc = x[0];\r
- double level = pars[0];\r
- double slope = pars[1];\r
- double slope2 = pars[2];\r
- double abscMin = pars[3];\r
- double abscMax = pars[4];\r
-\r
- if (absc >= abscMin && absc < abscMax) {\r
- TF1::RejectPoint();\r
- return 0;\r
- }\r
- \r
- double ordo = level + slope * absc +slope2 * absc*absc;\r
- return ordo;\r
-\r
-}\r
+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");
+}
--- /dev/null
+#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,"<N_{bin}> 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();
+}
--- /dev/null
+#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;j<NBinsArrays;j++){
+ signal[j]=0.0;
+ errSignal[j]=1.0;
+ } // <<< Not actually in use a the moment <<<
+
+ Stat_t mean[NBinsArrays], errMean[NBinsArrays]; // Mean and error
+ Stat_t sigma[NBinsArrays], errSigma[NBinsArrays]; // Gaussian width and error
+ Stat_t area[NBinsArrays], errArea[NBinsArrays]; // Peak area from fit and error
+ Stat_t integLow[NBinsArrays], integUp[NBinsArrays]; // Lower/upper limits for integration
+
+ // Float_t DefaultBinPtEdges[NDefault]={0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0};
+ Float_t BinPtEdges[NBinsArrays];
+
+ // Float_t BinPtEdges[NBins+1]={0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0};
+
+ Int_t BinLo, BinHi; //For pt bins
+ // **** Main loop over the FitControllers ****
+ TIter controlIter(ControlArray);
+ FitControl *controller;
+ while (controller=(FitControl*)controlIter.Next()) {
+ controller->CalcBinLimits(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;N<NBins-NFirst;N++){
+ NArray=N+1; // Arrays numbered from 1 not 0..
+ c1->cd(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;jj<NBinsArrays;jj++){
+// cout << "BinPtEdges " << jj << " = " << BinPtEdges[jj] << endl;
+// cout << "Mean " << jj << " = " << mean[jj] << endl;
+// cout << "Sigma " << jj << " = " << sigma[jj] << endl;
+// cout << "Signal + bkgd " << jj << " = " << sigBkgd[jj] << endl;
+// cout << "Background " << jj << " = " << bkgd[jj] << endl;
+// }
+ // Yields in each bin returned in a histogram
+ TH1F* hYields= new TH1F("Yield","Raw Particle Yield",NBinsArrays-1,BinPtEdges);
+ //hYields->Reset("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;
+
+}
debug=kTRUE
option="SAVE" #FIXME:set option
suffix=""
+fitFolder="LHC10h_000137161_p1_5plus"
+fitBin="00"
+partID=1
+task=no
+fit=no
+
give_help() {
This scripts runs the the physics selection and centrality on a specified run/dataset
Available options:
- Mode control, at least one of the following options should be used
+ * Run the task *
-r <mode> Run the task
Modes (compulsory):
0 local
(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 <suffix> Add extra suffix to files
Grid only options
-g <gridmode> Plugin Mode [default=$mode]
-p <recopass> Reconstruction pass [default=$pass]
-p <path> Data set path
-n <nev> Number of events
-w <workers> Number of workers [default=$workers]
+
+ * Fit the results *
+ -f <folder> Run the fitting macro in the subfolder of ./output
+ -b <bin> Centrality bin index [default=$fitBin]
+ -p <particleID> Fit particle defined by particleID [default=$partID]
+ 0=K0
+ 1=Lambda
+ 2=Anti-Lambda
+ 3=Lambda + Anti-Lambda,
+ 4=Csi
+ 6=Omega
+ -x <suffix> 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
;;
p)
pass=$OPTARG
+ partID=$OPTARG
+ ;;
+ x)
+ suffix=$OPTARG
;;
h)
give_help
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
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff