Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / macros / draw_separation.C
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
20 using 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
37 TF1* piFunc = 0;
38 TF1* kFunc  = 0;
39 TF1* pFunc = 0;
40 TF1* sigmaFunc = 0;
41
42 Double_t Sep(Double_t* xx, Double_t* par);
43
44 //____________________________________________________________________________
45 void 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 //______________________________________________________________________________
125 Double_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 }