generate and apply corrections fo PbPb analysis (Michael Knichel)
[u/mrichter/AliRoot.git] / PWG0 / dNdPt / macros / ApplyCorrections.C
CommitLineData
7c96ea29 1Double_t ApplyCorrections(const char* datafile, const char* datatask, const char* corrfile, const char* idstring ,const char* outfile, const char* gifdir = 0)
2//void ApplyCorrections()
3{
4
5// tmp setting
6// const char* datafile = "/lustre/alice/train/V005.PbPb/2010-11-28_0258.4230/mergedPeriods/data/PbPb/LHC10h.pass1/mknichel_dNdPtPbPb_gt_v0_c0.root";
7// const char* datatask = "mknichel_dNdPtPbPb_gt_v0_c0";
8// const char* corrfile = "/u/mknichel/alice/dNdPt/2010-11-28/corrMatr_gt_v0_c0.root";
9// const char* idstring = "gt_v0_c0";
10// const char* outfile = "/u/mknichel/alice/dNdPt/2010-11-28/finalSpectra_v0_c0.root";
11// const char* gifdir = "/u/mknichel/alice/dNdPt/2010-11-28";
12
13 // settings vor zVertex cut (event and track level)
14 Double_t zVert = 10.0;
15
16 // setting on eta cut (track level)
17 Double_t eta = 0.8;
18
19 //load required libraries
20 //load required libraries
21 gSystem->AddIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/ -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/PWG0 -I$ALICE_ROOT/PWG1 -I$ALICE_ROOT/PWG2 -I$ALICE_ROOT/PWG3 -I$ALICE_ROOT/PWG3/vertexingHF -I$ALICE_ROOT/PWG4 -I$ALICE_ROOT/CORRFW -I$ALICE_ROOT/TPC -I$ALICE_ROOT/TRD -I$ALICE_ROOT/PWG3/muon -I$ALICE_ROOT/JETAN -I$ALICE_ROOT/ANALYSIS/Tender");
22
23 gSystem->Load("libCore");
24 gSystem->Load("libPhysics");
25 gSystem->Load("libMinuit");
26 gSystem->Load("libGui");
27 gSystem->Load("libXMLParser");
28
29 gSystem->Load("libGeom");
30 gSystem->Load("libVMC");
31
32 gSystem->Load("libNet");
33
34 gSystem->Load("libSTEERBase");
35 gSystem->Load("libESD");
36 gSystem->Load("libCDB");
37 gSystem->Load("libRAWDatabase");
38 gSystem->Load("libRAWDatarec");
39 gSystem->Load("libANALYSIS");
40
41
42
43 gSystem->Load("libANALYSIS.so");
44 gSystem->Load("libANALYSISalice.so");
45 gSystem->Load("libTENDER.so");
46 gSystem->Load("libCORRFW.so");
47 gSystem->Load("libPWG0base.so");
48 gSystem->Load("libPWG0dep");
49 gSystem->Load("libPWG0selectors.so");
50
51
52 // make plots nicer
53 gROOT->SetStyle("Plain");
54 gStyle->SetPalette(1);
55
56 // array for all histograms to be saves
57 TObjArray* Hists = new TObjArray();
58
59 // open file with correction matrices
60 TFile *fcorr = TFile::Open(corrfile,"READ");
61 if (!fcorr) return -2;
62
63 // load data
64 TFile* fdata = TFile::Open(datafile,"READ");
65 if (!fdata) return -1;
66 TList* ldata = dynamic_cast<TList*>(fdata->Get(datatask));
67 if (!ldata) return -1;
68 AlidNdPtAnalysisPbPb *obj = dynamic_cast<AlidNdPtAnalysisPbPb*>(ldata->FindObject("dNdPtAnalysisPbPb"));
69 if (!obj) return -1;
70
71 //Event statistics
72 THnSparse *fRecEventHist2 = obj->GetRecEventHist2(); //reconstructed events
73 TH2D* h2RecEvent2All = fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2All");
74 fRecEventHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
75 TH2D* h2RecEvent2 = fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2");
76
77 Double_t ReconstructedEvents = h2RecEvent2->Integral();
78 Double_t ReconstructedEventsAll = h2RecEvent2All->Integral();
79
80 Hists->Add(h2RecEvent2);
81 Hists->Add(h2RecEvent2All);
82
83 printf("=== Number of events from DATA %lf ===\n",ReconstructedEvents);
84 printf("=== Number of events from DATA (before zVertex cut) %lf ===\n",ReconstructedEventsAll);
85
86 TH1D* h1ReconstructedEvents = new TH1D("h1ReconstructedEvents","h1ReconstructedEvents",1,0,1);
87 TH1D* h1ReconstructedEventsAll = new TH1D("h1ReconstructedEventsAll","h1ReconstructedEventsAll",1,0,1);
88
89 h1ReconstructedEvents->Fill(0,ReconstructedEvents);
90 h1ReconstructedEvents->SetEntries(ReconstructedEvents);
91 h1ReconstructedEventsAll->Fill(0,ReconstructedEventsAll);
92 h1ReconstructedEventsAll->SetEntries(ReconstructedEventsAll);
93
94 Hists->Add(h1ReconstructedEvents);
95 Hists->Add(h1ReconstructedEventsAll);
96
97 // retrieve tracks
98 THnSparse* fRecTrackHist2 = obj->GetRecTrackHist2(2); //after all cuts (2)
99 fRecTrackHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01); //zVertex
100 fRecTrackHist2->GetAxis(2)->SetRangeUser(-eta, eta-0.01); // eta
101
102 TH1D* h1RecTrackHist2_zv = fRecTrackHist2->Projection(0)->Clone("h1RecTrackHist2_zv");
103 TH1D* h1RecTrackHist2_pt = fRecTrackHist2->Projection(1)->Clone("h1RecTrackHist2_pt");
104 TH1D* h1RecTrackHist2_eta = fRecTrackHist2->Projection(2)->Clone("h1RecTrackHist2_eta");
105
106 Hists->Add(h1RecTrackHist2_zv);
107 Hists->Add(h1RecTrackHist2_pt);
108 Hists->Add(h1RecTrackHist2_eta);
109
110 // retrieve correction matrices for tracking efficiency (note different binning!)
111 TH1D* h1TrackCorr_pt = (TH1D*)fcorr->Get("h1TrackCorr_pt");
112 TH1D* h1TrackCorr_eta = (TH1D*)fcorr->Get("h1TrackCorr_eta");
113
114 // retrieve correction matrices for secondaries (note different binning!)
115 TH1D* h1SecCorr_pt = (TH1D*)fcorr->Get("h1SecCorr_pt");
116 TH1D* h1SecCorr_eta = (TH1D*)fcorr->Get("h1SecCorr_eta");
117
118 // create corrected spectra (as clone of raw data)
119 TH1D* h1Corrected_pt = h1RecTrackHist2_pt->Clone("h1Corrected_pt");
120 TH1D* h1Corrected_eta = h1RecTrackHist2_eta->Clone("h1Corrected_eta");
121
122 // secondaries correction for pt spectrum
123 for (int i=1; i <= h1Corrected_pt->GetNbinsX() ; i++) {
124 Double_t pt = h1Corrected_pt->GetBinCenter(i);
125 Double_t val = h1Corrected_pt->GetBinContent(i);
126 Double_t err = h1Corrected_pt->GetBinError(i);
127 if (pt >= 19) { pt = 19; } // above 20 GeV corr matr have low statistics
128 Double_t secCorr = h1SecCorr_pt->GetBinContent(h1SecCorr_pt->FindBin(pt));
129 Double_t secCorrErr = h1SecCorr_pt->GetBinError(h1SecCorr_pt->FindBin(pt));
130 Double_t cval = val*secCorr;
131 Double_t cerr = TMath::Sqrt(val*val*secCorrErr*secCorrErr + err*err*secCorr*secCorr);
132 h1Corrected_pt->SetBinContent(i,cval);
133 h1Corrected_pt->SetBinError(i,cerr);
134 }
135
136 // tracking efficiency correction pt spectrum
137 for (int i=1; i <= h1Corrected_pt->GetNbinsX() ;i++) {
138 Double_t pt = h1Corrected_pt->GetBinCenter(i);
139 Double_t val = h1Corrected_pt->GetBinContent(i);
140 Double_t err = h1Corrected_pt->GetBinError(i);
141 if (pt >= 19) { pt = 19; } // above 20 GeV corr matr have low statistics
142 Double_t effCorr = h1TrackCorr_pt->GetBinContent(h1TrackCorr_pt->FindBin(pt));
143 Double_t effCorrErr = h1TrackCorr_pt->GetBinError(h1TrackCorr_pt->FindBin(pt));
144
145 Double_t cval = val*effCorr;
146 Double_t cerr = TMath::Sqrt(val*val*effCorrErr*effCorrErr + err*err*effCorr*effCorr);
147 h1Corrected_pt->SetBinContent(i,cval);
148 h1Corrected_pt->SetBinError(i,cerr);
149 }
150
151 // for eta the correction is simpler because of same binning
152 h1Corrected_eta->Multiply(h1SecCorr_eta);
153 h1Corrected_eta->Multiply(h1TrackCorr_eta);
154
155 Hists->Add(h1Corrected_pt);
156 Hists->Add(h1Corrected_eta);
157
158 // create final spectra (as clone of corrected data)
159 TH1D* dNdPt = h1Corrected_pt->Clone("dNdPt");
160 TH1D* dNdEta = h1Corrected_eta->Clone("dNdEta");
161
162 // also uncorrected spectra (as clone of raw data)
163 TH1D* dNdPt_raw = h1RecTrackHist2_pt->Clone("dNdPt_raw");
164 TH1D* dNdEta_raw = h1RecTrackHist2_eta->Clone("dNdEta_raw");
165
166
167//TF1 *fperi = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,40.);
168//TF1 *fcent = new TF1("cent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,40.);
169TString id = TString(idstring);
170if ( id.Contains("c0") || id.Contains("c5") ) {
171 TF1 *fun = new TF1("cent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,40.);
172 cout << "++++++++++++++++++++++++++++++++++++++" << endl;
173 cout << "CENTRAL pt-resolution correction used!" << endl;
174 cout << "++++++++++++++++++++++++++++++++++++++" << endl;
175} else {
176TF1 *fun = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,40.);
177 cout << "+++++++++++++++++++++++++++++++++++++++++" << endl;
178 cout << "PERIPHERAL pt-resolution correction used!" << endl;
179 cout << "+++++++++++++++++++++++++++++++++++++++++" << endl;
180}
181
182 // normalization and finalization
183 // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
184 for (int i=1; i <= dNdPt->GetNbinsX() ;i++) {
185 Double_t pt = dNdPt->GetBinCenter(i);
186 Double_t width = dNdPt->GetBinWidth(i);
187 Double_t val = dNdPt->GetBinContent(i);
188 Double_t err = dNdPt->GetBinError(i);
189 Double_t corrPtResol = 1.0;
190 corrPtResol = fun->Eval(pt);
191 if (pt < 5.) { corrPtResol = 1.0; }
192 Double_t cval = (val * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
193 Double_t cerr = (err * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
194 dNdPt->SetBinContent(i,cval);
195 dNdPt->SetBinError(i,cerr);
196 }
197 // for dndeta again simpler
198 dNdEta->Scale(1,"width");
199 dNdEta->Scale(1./ReconstructedEvents);
200
201 // normalization and finalization
202 // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
203 for (int i=1; i <= dNdPt_raw->GetNbinsX() ;i++) {
204 Double_t pt = dNdPt_raw->GetBinCenter(i);
205 Double_t width = dNdPt_raw->GetBinWidth(i);
206 Double_t val = dNdPt_raw->GetBinContent(i);
207 Double_t err = dNdPt_raw->GetBinError(i);
208 Double_t cval = val/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
209 Double_t cerr = err/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
210 dNdPt_raw->SetBinContent(i,cval);
211 dNdPt_raw->SetBinError(i,cerr);
212 }
213 // for dndeta again simpler
214 dNdEta_raw->Scale(1,"width");
215 dNdEta_raw->Scale(1./ReconstructedEvents);
216
217 dNdEta->SetMarkerStyle(21);
218 dNdPt->SetMarkerStyle(21);
219 dNdPt->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
220 dNdEta->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
221
222 dNdPt_raw->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
223 dNdEta_raw->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
224
225 Hists->Add(dNdEta_raw);
226 Hists->Add(dNdPt_raw);
227 Hists->Add(dNdEta);
228 Hists->Add(dNdPt);
229
230
231 // plot pictures and save to gifdir
232 for (i=0; i < Hists->LastIndex(); i++) {
233 TCanvas* ctmp = PlotHist(Hists->At(i),idstring);
234 if (gifdir && ctmp) {
235 TString gif(gifdir);
236 gif += '/';
237 gif += ctmp->GetName();
238 gif += ".gif";
239 ctmp->SaveAs(gif.Data(),"gif");
240 delete ctmp;
241 }
242 }
243
244 // save all correction matrices and control histograms to file
245 if (!outfile) { return; }
246 TFile *out = TFile::Open(outfile,"RECREATE");
247 Hists->Write();
248 out->Close();
249
250 return ReconstructedEvents;
251
252}
253
254//___________________________________________________________________________
255TCanvas* PlotHist(TObject* hobj, const char* label=0)
256{
257 TH1* h = dynamic_cast<TH1*>(hobj);
258 if (!h) return 0;
259 if (h->GetDimension() > 2) return 0;
260 h->SetStats(0);
261 if ( TString(h->GetName()).Contains("Events")) { h->SetStats(1); }
262 TString t(label);
263 if (label) t += "_";
264 t += h->GetName();
265 h->SetTitle(t.Data());
266 TCanvas* c = new TCanvas(t.Data(),t.Data());
267 if (h->GetDimension() >= 1) {
268 TString xlabel(h->GetXaxis()->GetTitle());
269 if (xlabel.Contains("Pt")) { c->SetLogx(); c->SetLogy(); h->GetXaxis()->SetRangeUser(0.1 , 50.); }
270 if (xlabel.Contains("p_{T}")) { c->SetLogx(); c->SetLogy(); h->GetXaxis()->SetRangeUser(0.1 , 50.); }
271 }
272 if (h->GetDimension() == 2) {
273 TString ylabel(h->GetYaxis()->GetTitle());
274 if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 50.); }
275 if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 50.); }
276 h->Draw("COLZ");
277 }
278 if (h->GetDimension() == 1) {
279 h->Draw();
280 }
281 return c;
282
283}