Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / macros / draw_separation.C
CommitLineData
4ebdd20e 1#include <TFile.h>
2#include <TList.h>
3#include <TTree.h>
4#include <TChain.h>
5#include <TClonesArray.h>
6#include <TObjString.h>
7#include <TString.h>
8#include <TROOT.h>
9#include <TLegend.h>
10#include <TStyle.h>
11#include <TCanvas.h>
12
13#include "my_tools.C"
14#include "my_functions.C"
15
16#include <iostream>
17#include <fstream>
18#include <string>
19
20using namespace std;
21
22/*
23 To run calibrations:
24 ====================
25
26 Use root:
27
28 .L libMyDeDxAnalysis.so
29 .L my_functions.C+
30 .L my_tools.C+
31 .L draw_separation.C+
32
33 DrawSeparation("fitparameters/lhc10h_aod_all.root", 0, 50)
34 */
35
36
37TF1* piFunc = 0;
38TF1* kFunc = 0;
39TF1* pFunc = 0;
40TF1* sigmaFunc = 0;
41
42Double_t Sep(Double_t* xx, Double_t* par);
43
44//____________________________________________________________________________
45void DrawSeparation(const Char_t* fitFileName, Double_t pLow, Double_t pHigh)
46{
47 gStyle->SetOptStat(0);
48
49 TFile* fitFile = FindFileFresh(fitFileName);
50 if(!fitFile)
51 return;
52 DeDxFitInfo* fitPar = (DeDxFitInfo*)fitFile->Get("fitInfo");
53 fitPar->Print();
54
55 fixMIP = fitPar->MIP;
56 fixPlateau = fitPar->plateau;
57
58 Double_t dedxPar[6] = {0, 0, 0, 0, 0, 0};
59 Double_t sigmaPar[6] = {0, 0, 0, 0, 0, 0};
60
61 dedxPar[0] = fitPar->optionDeDx;
62 for(Int_t i = 0; i < fitPar->nDeDxPar; i++) {
63 dedxPar[i+1] = fitPar->parDeDx[i];
64 }
65
66 sigmaPar[0] = fitPar->optionSigma;
67 for(Int_t i = 0; i < fitPar->nSigmaPar; i++) {
68 sigmaPar[i+1] = fitPar->parSigma[i];
69 }
70
71 piFunc = new TF1("piFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
72 piFunc->SetParameters(dedxPar);
73
74 kFunc = new TF1("kFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
75 kFunc->SetParameters(dedxPar);
76 kFunc->SetParameter(0, kFunc->GetParameter(0)+10);
77
78 pFunc = new TF1("pFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
79 pFunc->SetParameters(dedxPar);
80 pFunc->SetParameter(0, pFunc->GetParameter(0)+20);
81
82 sigmaFunc = new TF1("sigmaFunc", SigmaFunc, 0, 100, fitPar->nSigmaPar+1);
83 sigmaFunc->SetParameters(sigmaPar);
84
85 TCanvas* c1 = new TCanvas("c1", "c1");
86
87 TH1F* hist = new TH1F("hist", "Separation in pp vs p; p [GeV/c]; Separation",
88 100, 0, pHigh);
89 hist->SetMinimum(0.0);
90 hist->SetMaximum(6.0);
91 hist->Draw();
92
93 TLegend* legend = new TLegend(0.74, 0.64, 0.89, 0.89);
94 legend->SetBorderSize(0);
95 legend->SetFillColor(0);
96
97 TF1* pipFunc = new TF1("pipFunc", Sep,
98 pLow, pHigh, 1);
99 pipFunc->SetParameter(0, 0);
100 pipFunc->SetLineColor(2);
101 pipFunc->Draw("SAME");
102
103 TF1* pikFunc = new TF1("pikFunc", Sep,
104 pLow, pHigh, 1);
105 pikFunc->SetParameter(0, 1);
106 pikFunc->SetLineColor(3);
107 pikFunc->Draw("SAME");
108
109 TF1* kpFunc = new TF1("kpFunc", Sep,
110 pLow, pHigh, 1);
111 kpFunc->SetParameter(0, 2);
112 kpFunc->SetLineColor(4);
113 kpFunc->Draw("SAME");
114
115 legend->AddEntry(pipFunc, "#pi-p", "L");
116 legend->AddEntry(pikFunc, "#pi-K", "L");
117 legend->AddEntry(kpFunc, "K-p", "L");
118 legend->Draw();
119 gROOT->ProcessLine(".x drawText.C");
120 c1->SaveAs("separation.gif");
121 c1->SaveAs("separation.pdf");
122}
123
124//______________________________________________________________________________
125Double_t Sep(Double_t* xx, Double_t* par)
126{
127 //
128 // Could speed up fit by forcing it to use <p>. In that way the parameters
129 // could be amde statis cand only changed when going to a new p bin
130 //
131 Double_t p = xx[0];
132
133 Int_t option = Int_t(par[0]);
134
135 TF1* f1 = 0;
136 TF1* f2 = 0;
137 switch (option) {
138
139 case 0: // pi - p
140 f1 = piFunc;
141 f2 = pFunc;
142 break;
143 case 1: // pi - k
144 f1 = piFunc;
145 f2 = kFunc;
146 break;
147 case 2: // k - p
148 f1 = kFunc;
149 f2 = pFunc;
150 break;
151 default:
152 cout << "Error in Sep: option " << option << " not supported!!!!!" << endl;
153 return 0;
154 break;
155 }
156
157 Double_t dedx1 = f1->Eval(p);
158 Double_t dedx2 = f2->Eval(p);
159 Double_t sigma1 = sigmaFunc->Eval(dedx1);
160 Double_t sigma2 = sigmaFunc->Eval(dedx2);
161
162 return (dedx1-dedx2)/TMath::Sqrt(sigma1*sigma2);
163}