]>
Commit | Line | Data |
---|---|---|
943dc129 | 1 | #include "TTree.h" |
2 | #include "TFile.h" | |
3 | #include <TList.h> | |
4 | #include <Rtypes.h> | |
5 | #include "TCanvas.h" | |
6 | #include <TH2I.h> | |
7 | #include <TH2F.h> | |
8 | #include <TProfile.h> | |
9 | #include <TFitResultPtr.h> | |
10 | #include <TFitResult.h> | |
11 | #include "TF1.h" | |
12 | #include <iostream> | |
13 | #include <AliAnalysisEtTrackMatchCorrections.h> | |
14 | #include <AliAnalysisEtRecEffCorrection.h> | |
15 | ||
16 | TCanvas *c1 = 0; | |
17 | ||
18 | TF1 * generateRecEffFunction(Double_t p0, Double_t p1); | |
19 | ||
20 | int calculateCorrections(TString filename="Et.ESD.simPbPb.PHOS.root", Double_t p0 = 0.366, Double_t p1 = 0.0) | |
21 | { | |
22 | c1 = new TCanvas; | |
23 | ||
24 | TFile *f = TFile::Open(filename, "READ"); | |
25 | ||
26 | TList *l = (TList*)f->Get("out1"); | |
27 | ||
28 | TTree *primaryTree = (TTree*)l->FindObject("fPrimaryTreePhosMC"); | |
29 | TTree *recTree = (TTree*)l->FindObject("fEventSummaryTreePhosRec"); | |
30 | TTree *mcTree = (TTree*)l->FindObject("fEventSummaryTreePhosMC"); | |
31 | std::cout << primaryTree << " " << recTree << " " << mcTree << std::endl; | |
32 | ||
33 | Int_t clusterMult = 0; | |
34 | Int_t nChargedNotRemoved = 0; | |
35 | Int_t nNeutralNotRemoved = 0; | |
36 | Int_t nGammaRemoved = 0; | |
37 | Int_t nSecNotRemoved = 0; | |
38 | ||
39 | Double_t emEtRec = 0; | |
40 | Double_t emEtMc = 0; | |
41 | ||
42 | recTree->SetBranchAddress("fNeutralMultiplicity", &clusterMult); | |
43 | mcTree->SetBranchAddress("fChargedNotRemoved", &nChargedNotRemoved); | |
44 | mcTree->SetBranchAddress("fNeutralNotRemoved", &nNeutralNotRemoved); | |
45 | mcTree->SetBranchAddress("fGammaRemoved", &nGammaRemoved); | |
46 | mcTree->SetBranchAddress("fSecondaryNotRemoved", &nSecNotRemoved); | |
47 | ||
48 | TH2I *hChargedVsClusterMult = new TH2I("hChVsMult", "hChVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5); | |
49 | TH2I *hNeutralVsClusterMult = new TH2I("hNeutVsMult", "hNeutVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5); | |
50 | TH2I *hGammaVsClusterMult = new TH2I("hGammaVsMult", "hGammaVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5); | |
51 | TH2I *hSecVsClusterMult = new TH2I("hSecVsMult", "hSecVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5); | |
52 | ||
53 | Int_t nEvents = mcTree->GetEntriesFast(); | |
54 | for(Int_t i = 0; i < nEvents; i++) | |
55 | { | |
56 | mcTree->GetEvent(i); | |
57 | recTree->GetEvent(i); | |
58 | hChargedVsClusterMult->Fill(clusterMult, nChargedNotRemoved); | |
59 | hNeutralVsClusterMult->Fill(clusterMult, nNeutralNotRemoved); | |
60 | hGammaVsClusterMult->Fill(clusterMult, nGammaRemoved); | |
61 | hSecVsClusterMult->Fill(clusterMult, nSecNotRemoved); | |
62 | } | |
63 | ||
64 | c1->Divide(2,2); | |
65 | c1->cd(1); | |
66 | TString title = "Charged particles not removed"; | |
67 | TString xtitle = "Cluster Multiplicity"; | |
68 | TString ytitle = "N_{ch}"; | |
69 | hChargedVsClusterMult->SetTitle(title); | |
70 | hChargedVsClusterMult->GetYaxis()->SetTitle(ytitle); | |
71 | hChargedVsClusterMult->GetXaxis()->SetTitle(xtitle); | |
72 | hChargedVsClusterMult->SetStats(0); | |
73 | hChargedVsClusterMult->Draw(); | |
74 | TProfile *chProf = hChargedVsClusterMult->ProfileX(); | |
75 | chProf->SetStats(0); | |
76 | chProf->Draw("SAME"); | |
77 | //TF1 fitcharged("fitcharged","([0] + [1]*x)*(0.48/([2] + [3]*[2]))", 0, 100); | |
78 | TF1 fitcharged("fitcharged","pol1", 0, 100); | |
79 | // fitcharged.FixParameter(2, p0); | |
80 | // fitcharged.FixParameter(3, p1); | |
81 | TFitResultPtr chRes = chProf->Fit(&fitcharged,"S"); | |
82 | TArrayD ch; | |
83 | if(!chRes) | |
84 | { | |
85 | ch = TArrayD(chRes->NPar(),chRes->GetParams()); | |
86 | } | |
87 | else | |
88 | { | |
89 | std::cout << "Could not extract charged contribution params" << std::endl; | |
90 | } | |
91 | c1->cd(2); | |
92 | title = "Neutral particles not removed"; | |
93 | ytitle = "N_{neutral}"; | |
94 | hNeutralVsClusterMult->SetTitle(title); | |
95 | hNeutralVsClusterMult->GetXaxis()->SetTitle(xtitle); | |
96 | hNeutralVsClusterMult->GetYaxis()->SetTitle(ytitle); | |
97 | hNeutralVsClusterMult->SetStats(0); | |
98 | hNeutralVsClusterMult->Draw(); | |
99 | TProfile *neuProf = hNeutralVsClusterMult->ProfileX(); | |
100 | neuProf->SetStats(0); | |
101 | neuProf->Draw("SAME"); | |
102 | TF1 fitneutral("fitneutral","pol1", 0, 100); | |
103 | TFitResultPtr neuRes = neuProf->Fit(&fitneutral,"S"); | |
104 | TArrayD neu; | |
105 | if(!neuRes) | |
106 | { | |
107 | neu = TArrayD(neuRes->NPar(),neuRes->GetParams()); | |
108 | } | |
109 | else | |
110 | { | |
111 | std::cout << "Could not extract charged contribution params" << std::endl; | |
112 | } | |
113 | c1->cd(3); | |
114 | ||
115 | title = "Gammas removed"; | |
116 | ytitle = "N_{#gamma}"; | |
117 | hGammaVsClusterMult->SetTitle(title); | |
118 | hGammaVsClusterMult->GetYaxis()->SetTitle(ytitle); | |
119 | hGammaVsClusterMult->GetXaxis()->SetTitle(xtitle); | |
120 | hGammaVsClusterMult->SetStats(0); | |
121 | hGammaVsClusterMult->Draw(); | |
122 | TProfile *gamProf = hGammaVsClusterMult->ProfileX(); | |
123 | gamProf->SetStats(0); | |
124 | gamProf->Draw("SAME"); | |
125 | TF1 fitgamma("fitgamma","pol1", 0, 100); | |
126 | TFitResultPtr gammaRes = gamProf->Fit(&fitgamma,"S"); | |
127 | TArrayD gamma; | |
128 | if(!gammaRes) | |
129 | { | |
130 | gamma = TArrayD(gammaRes->NPar(),gammaRes->GetParams()); | |
131 | } | |
132 | else | |
133 | { | |
134 | std::cout << "Could not extract charged contribution params" << std::endl; | |
135 | } | |
136 | ||
137 | c1->cd(4); | |
138 | ||
139 | title = "Secondaries not removed"; | |
140 | ytitle = "N_{sec}"; | |
141 | hSecVsClusterMult->SetTitle(title); | |
142 | hSecVsClusterMult->GetXaxis()->SetTitle(xtitle); | |
143 | hSecVsClusterMult->GetYaxis()->SetTitle(ytitle); | |
144 | hSecVsClusterMult->SetStats(0); | |
145 | hSecVsClusterMult->Draw(); | |
146 | TProfile *secProf = hSecVsClusterMult->ProfileX(); | |
147 | secProf->SetStats(0); | |
148 | secProf->Draw("SAME"); | |
149 | TF1 fitsecondary("fitsecondary","pol1", 0, 100); | |
150 | TFitResultPtr secRes = secProf->Fit(&fitsecondary,"S"); | |
151 | TArrayD sec; | |
152 | if(!secRes) | |
153 | { | |
154 | sec = TArrayD(secRes->NPar(),secRes->GetParams()); | |
155 | } | |
156 | else | |
157 | { | |
158 | std::cout << "Could not extract charged contribution params" << std::endl; | |
159 | } | |
160 | Double_t meanCharged = 0.48/(p0 + p1*0.48); | |
161 | Double_t meanNeutral = 0.53/(p0 + p1*0.53); | |
162 | Double_t meanGamma = 0.51/(p0 + p1*0.51); | |
163 | Double_t meanSecondary = meanGamma; | |
164 | ||
165 | AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections("TmCorrectionsPhos", fitcharged, fitneutral, fitgamma, fitsecondary, | |
166 | meanCharged, meanNeutral, meanGamma, meanSecondary ); | |
167 | ||
168 | TF1 *func = generateRecEffFunction(p0, p1); | |
169 | AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection("ReCorrectionsPhos", *func, 1000); | |
170 | TFile *outfile = TFile::Open("calocorrections.root","RECREATE"); | |
171 | cor->Write(); | |
172 | recor->Write(); | |
173 | return 0; | |
174 | } | |
175 | ||
176 | ||
177 | TF1* generateRecEffFunction(Double_t p0, Double_t p1) | |
178 | { | |
179 | TF1 *f = new TF1("receff", "x/([0] + x*[1])", 0, 200); | |
180 | ||
181 | Double_t params[2] = {p0, p1}; | |
182 | f->SetParameters(params); | |
183 | return f; | |
184 | } | |
185 | ||
186 |