Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / final_spectra / makeSpectra.C
CommitLineData
4ebdd20e 1#include <TFile.h>
2#include <TH1.h>
3#include <TF1.h>
4#include <TCanvas.h>
5#include <TList.h>
6#include <TLegend.h>
7#include <TStyle.h>
8#include <TMath.h>
9#include <TSystem.h>
10#include <TLatex.h>
11#include <TASImage.h>
12#include <TLine.h>
13#include <TBox.h>
14
15#include "my_tools.C"
16
17
18/*
19 Info:
20 Multiplies the fraction with the charged spectrum.
21 Applies a rapidity correction in case of K and p.
22
23 To run:
24
25 gSystem->AddIncludePath("-I../macros")
26 gSystem->AddIncludePath("-I../lib")
27 gROOT->SetMacroPath(".:../macros")
28 .L my_tools.C+
29 .L makeSpectra.C+
30
31 MakeSpectra("../corrected_fraction/fraction_7tev_b_all.root", "NOPID_pp_7TeV.root", "final_spectra_7tev_b.root");
32
33
34 */
35
36Double_t rap_correction(Double_t* x, Double_t* par);
37
38void MakeSpectra(const Char_t* fileName,
39 const Char_t* fileNameCharged,
40 const Char_t* outFileName,
41 const Char_t* endName="PP")
42{
43 TFile* fileData = FindFileFresh(fileName);
44
45 if(!fileData)
46 return;
47
48 TFile* fileDataCharged = FindFileFresh(fileNameCharged);
49
50 if(!fileDataCharged)
51 return;
52
53 const Int_t nPid = 3;
54 const Char_t* pidNames[nPid] = {"Pion", "Kaon", "Proton"};
55 const Char_t* titleNames[nPid] = {"#pi^{+} + #pi^{-}", "K^{+} + K^{-}", "p + #bar{p}"};
56 const Double_t mass[nPid] = {0.140, 0.494, 0.938};
57
58 TH1D* hPtSpectrum[nPid] = {0, 0, 0};
59 TH1D* hPtSpectrumSyst[nPid] = {0, 0, 0};
60
61 TH1D* hPtCharged = (TH1D*)fileDataCharged->Get(Form("hDnDptCharged_%s", endName));
62 TH1D* hPtChargedSyst = (TH1D*)fileDataCharged->Get(Form("hDnDptChargedSyst_%s", endName));
63
64 TF1* fRap = new TF1("fRap", rap_correction, 0.0, 50.0, 2);
65
66 for(Int_t pid = 0; pid < nPid; pid++) {
67
68 // Pions
69 TH1D* hFraction = (TH1D*)fileData->Get(Form("h%sFraction_%s", pidNames[pid], endName));
70 TH1D* hFractionSyst = (TH1D*)fileData->Get(Form("h%sFractionSyst_%s", pidNames[pid], endName));
71
72 hPtSpectrum[pid] = (TH1D*)hFraction->Clone(Form("h%sSpectrum_%s", pidNames[pid], endName));
73 hPtSpectrum[pid]->GetYaxis()->SetTitle(Form("1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{%s})/(dy dp_{T}) (GeV/c)^{-2}", titleNames[pid]));
74 hPtSpectrum[pid]->Multiply(hPtCharged);
75
76 hPtSpectrumSyst[pid] = (TH1D*)hFractionSyst->Clone(Form("h%sSpectrumSyst_%s", pidNames[pid], endName));
77 hPtSpectrumSyst[pid]->SetMarkerStyle(1);
78 hPtSpectrum[pid]->GetYaxis()->SetTitle(Form("1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{%s})/(dy dp_{T}) (GeV/c)^{-2}", titleNames[pid]));
79 hPtSpectrumSyst[pid]->Multiply(hPtChargedSyst);
80 hPtSpectrumSyst[pid]->SetFillStyle(1001);
81 hPtSpectrumSyst[pid]->SetFillColor(kGray);
82
83 //
84 // Rapidity correction
85 //
86 fRap->SetParameters(0.8, mass[pid]);
87 hPtSpectrum[pid]->Divide(fRap);
88 hPtSpectrumSyst[pid]->Divide(fRap);
89 }
90
91 TFile* fileOut = new TFile(outFileName, "RECREATE");
92 for(Int_t pid = 0; pid < nPid; pid++) {
93 hPtSpectrumSyst[pid]->Write();
94 hPtSpectrum[pid]->Write();
95 }
96 fileOut->Close();
97}
98
99//______________________________________________________________________________
100Double_t rap_correction(Double_t* x, Double_t* par)
101{
102 Double_t pt = x[0];
103
104 Double_t eta = par[0];
105 Double_t mass = par[1];
106
107 const Double_t mt = TMath::Sqrt(pt*pt + mass*mass);
108
109 const Double_t rap = TMath::ASinH(pt/mt*TMath::SinH(eta));
110
111 return rap/eta;
112}