]>
Commit | Line | Data |
---|---|---|
29e08f93 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | ||
18 | // Macro to compare/plot the Ntuple stored in MUONefficiency.root | |
53399039 | 19 | // comparison is done between the generated and reconstructed resonance |
20 | // Default is Upsilon but Jpsi can be studied if the ResType argument is changed | |
21 | // This allows to determine several important quantities | |
22 | // reconstruction efficiency (versus pt,y), invariant mass peak variations (vs pt,y) | |
29e08f93 | 23 | // reconstructed tracks and trigger tracks matching efficiency |
24 | ||
25 | // Christophe Suire, IPN Orsay | |
26 | ||
53399039 | 27 | void MUONplotefficiency(Int_t ResType = 553, Int_t fittype = 1){ |
29e08f93 | 28 | |
29 | Int_t SAVE=1; | |
30 | Int_t WRITE=1; | |
53399039 | 31 | |
32 | Double_t MUON_MASS = 0.105658369; | |
33 | Double_t UPSILON_MASS = 9.4603 ; | |
34 | Double_t JPSI_MASS = 3.096916; | |
35 | Double_t PI = 3.14159265358979312; | |
36 | Double_t RESONANCE_MASS; | |
37 | ||
38 | if (ResType==553) | |
39 | RESONANCE_MASS = UPSILON_MASS; | |
40 | if (ResType==443) | |
41 | RESONANCE_MASS = JPSI_MASS ; | |
42 | ||
43 | // 553 is the pdg code for Upsilon | |
44 | // 443 is the pdg code for Jpsi | |
29e08f93 | 45 | |
46 | Int_t fitfunc = fittype; | |
47 | // 0 gaussian fit +/- 450 MeV/c2 | |
48 | // 1 gaussian fit +/- 250 MeV/c2 | |
53399039 | 49 | // the following are not implemented in this version |
29e08f93 | 50 | // 2 approximated landau fit reverted |
51 | // 3 landau fit reverted convoluated with a gaussian | |
52 | // 4 reverted landau fitlanUpsilon | |
53 | ||
53399039 | 54 | |
29e08f93 | 55 | Float_t gaussianFitWidth = 0.450 ; // in GeV/c2 |
56 | if (fitfunc==1) gaussianFitWidth = 0.250 ; | |
53399039 | 57 | |
29e08f93 | 58 | // fit range : very important for LandauGauss Fit |
53399039 | 59 | Float_t FitLow = 0.0 ; Float_t FitHigh = 100. ; |
60 | if (ResType==443) FitLow = 2.2 ; FitHigh = 3.9 ; | |
61 | if (ResType==553) FitLow = 8.5 ; FitHigh = 10.2 ; | |
62 | ||
29e08f93 | 63 | |
53399039 | 64 | //gStyle->Reset(); |
65 | // two ways to set the stat box for all histos | |
29e08f93 | 66 | gStyle->SetOptStat(1110); |
67 | ||
68 | // place the stat box | |
69 | gStyle->SetStatX(0.5); | |
70 | gStyle->SetStatY(0.550); | |
71 | //gStyle->SetStatW(0.2); | |
72 | gStyle->SetStatH(0.2); | |
29e08f93 | 73 | |
74 | // choose histos with stat | |
75 | // gStyle->SetOptStat(1111); | |
76 | //hist->SetStats(kTRUE); or hist->SetStats(kFALSE); | |
77 | ||
78 | //gStyle->SetOptFit(1111); | |
79 | gStyle->SetOptFit(1); | |
80 | //gStyle->SetOptFit(0); | |
81 | ||
82 | gStyle->SetOptTitle(0); | |
83 | gStyle->SetPalette(1,0); | |
84 | //gStyle->SetOptLogy(1); | |
85 | //gStyle->SetOptLogx(0); | |
86 | ||
87 | gStyle->SetFrameFillColor(10); | |
88 | gStyle->SetFillColor(10); | |
89 | gStyle->SetCanvasColor(10); | |
90 | ||
91 | ||
92 | TFile *effFile = new TFile("MUONefficiency.root","READ"); | |
29e08f93 | 93 | TNtuple *Ktuple = (TNtuple*)effFile->Get("Ktuple"); |
53399039 | 94 | //("Ktuple","Kinematics NTuple","ev:npart:id:idmo:idgdmo:p:pt:y:theta:pseudorap:vx:vy:vz"); |
29e08f93 | 95 | TNtuple *ESDtuple = (TNtuple*)effFile->Get("ESDtuple"); |
53399039 | 96 | //("ESDtuple","Reconstructed Mu+Mu- pairs","ev:pt:y:theta:minv:pt1:y1:theta1:q1:pt2:y2:theta2:q2"); |
97 | ||
98 | ||
29e08f93 | 99 | /*********************************/ |
100 | // Histograms limits and binning | |
101 | Float_t ptmin = 0.0; Float_t ptmax = 20.0; Int_t ptbins = 10; | |
102 | Float_t ymin = -4.5; Float_t ymax = -2.; Int_t ybins = 10; | |
103 | Float_t thetamin = 165.; Float_t thetamax = 180.; Int_t thetabins = 100; | |
104 | ||
105 | Float_t etacutmin = -4.04813 ; | |
106 | Float_t etacutmax = -2.54209 ; | |
53399039 | 107 | |
108 | Float_t invMassMin = .0 ; Float_t invMassMax = 15.0 ; Int_t invMassBins = 100 ; | |
109 | Float_t ptmin = 0.0; Float_t ptmax = 20.0; Int_t ptbins = 10; | |
110 | if (ResType==443){ | |
111 | invMassMin = 1.0 ; invMassMax = 4.5 ; | |
112 | ptmin = 0.0; ptmax = 10.0; | |
113 | } | |
114 | if (ResType==553){ | |
115 | invMassMin = 7.0 ; invMassMax = 11.0 ; | |
116 | ptmin = 0.0; ptmax = 20.0; | |
117 | } | |
29e08f93 | 118 | |
119 | /*********************************/ | |
120 | // Values used to define the acceptance | |
121 | // | |
122 | Float_t thetacutmin = 171.0; | |
123 | Float_t thetacutmax = 178.; | |
124 | ||
125 | Float_t ptcutmin = 0.; | |
126 | Float_t ptcutmax = 20.; | |
127 | ||
53399039 | 128 | if (ResType==443){ |
129 | ptcutmin = 0.; ptcutmax = 10.; | |
130 | } | |
131 | ||
132 | if (ResType==553){ | |
133 | ptcutmin = 0.; ptcutmax = 20.; | |
134 | } | |
135 | ||
136 | ||
29e08f93 | 137 | Float_t masssigma = 0.1; // 100 MeV/c2 is a correct estimation |
138 | Float_t masscutmin = 0; Float_t masscutmax = 0 ; | |
139 | if (masssigma){ | |
53399039 | 140 | masscutmin = RESONANCE_MASS - 3.0*masssigma ; |
141 | masscutmax = RESONANCE_MASS + 3.0*masssigma ; | |
29e08f93 | 142 | } |
143 | else { | |
53399039 | 144 | masscutmin = RESONANCE_MASS - 1.0 ; |
145 | masscutmax = RESONANCE_MASS + 1.0 ; | |
29e08f93 | 146 | } |
147 | ||
148 | ||
29e08f93 | 149 | // here no cut on theta is necesary since during the simulation |
150 | // gener->SetChildThetaRange(171.0,178.0); | |
53399039 | 151 | // thus the resonance are generated in 4 PI but only the ones for which the decay muons |
29e08f93 | 152 | // are in the dimuon arm acceptance |
53399039 | 153 | // the pt cut is also "critical", in the generation part, resonance are generated within the range |
154 | // (ptcutmin,ptcutmax). During the reconstruction, the resonance could have a pt greater than ptcutmax !! | |
155 | // Should these resonances be cut or not ? | |
29e08f93 | 156 | // probably not but they will be for now (~ 1/2000) |
157 | ||
158 | // Another acceptance cut to add is a range for the recontructed invariant mass, it is | |
159 | // obvious that an upsilon reconstructed with a mass of 5 GeV/c2 is not correct. Thus | |
53399039 | 160 | Char_t ResonanceAccCutMC[200]; |
161 | sprintf(ResonanceAccCutMC,"pt>= %.2f && pt<= %.2f",ptcutmin,ptcutmax); | |
162 | ||
163 | Char_t ResonanceAccCutESD[200]; | |
164 | sprintf(ResonanceAccCutESD,"pt >= %.2f && pt<= %.2f && minv>= %.2f && minv<= %.2f",ptcutmin,ptcutmax,masscutmin,masscutmax); | |
165 | ||
29e08f93 | 166 | |
29e08f93 | 167 | |
168 | /*********************************/ | |
169 | // Cut conditions (Id,Background...) | |
170 | ||
53399039 | 171 | Char_t IdcutResonanceMC[100]; |
172 | Char_t IdcutMuminusMC[100]; | |
173 | Char_t IdcutMuplusMC[100]; | |
174 | ||
175 | if (ResType==553){ | |
176 | sprintf(IdcutResonanceMC,"id==553"); | |
177 | sprintf(IdcutMuminusMC,"id==13 && idmo==553"); | |
178 | sprintf(IdcutMuplusMC,"id==-13 && idmo==553"); | |
179 | } | |
180 | if (ResType==443){ | |
181 | sprintf(IdcutResonanceMC,"id==443"); | |
182 | sprintf(IdcutMuminusMC,"id==13 && idmo==443"); | |
183 | sprintf(IdcutMuplusMC,"id==-13 && idmo==443"); | |
184 | } | |
185 | ||
29e08f93 | 186 | //means no cuts since we don't have the trackid propagated yet pt>0 |
187 | // now we have it 10/05/2005 but it's still not being used | |
188 | ||
189 | ||
53399039 | 190 | // Background is calculated in the MUONefficiency.C macro |
191 | // Background calculation is meaningful only when Resonance have been generated with realistic background | |
192 | // when it's a pure Resonance file, the background lies artificially at the Resonance mass | |
29e08f93 | 193 | |
194 | //no realistic background | |
53399039 | 195 | Int_t REALISTIC_BACKGROUND = 0 ; |
29e08f93 | 196 | |
53399039 | 197 | // to take background into account, i.e. substraction, if necessary |
198 | // this is not ready yet | |
199 | Char_t BckgdCutResonanceESD[100]; | |
200 | sprintf(BckgdCutResonanceESD,"pt>0 && minv>%f && minv<%f && pt1>0 && pt2>0",masscutmin,masscutmax); | |
201 | ||
202 | // if realistic background | |
29e08f93 | 203 | // same cut + substract the background from hInvMassBg |
204 | ||
29e08f93 | 205 | |
206 | ||
53399039 | 207 | /*******************************************************************************************************************/ |
29e08f93 | 208 | Char_t txt[50]; |
209 | TLatex *tex; | |
210 | TLegend *leg; | |
29e08f93 | 211 | |
212 | /*******************************/ | |
213 | /* Monte Carlo Part */ | |
214 | /*******************************/ | |
215 | ||
216 | //---------------------------- | |
217 | // Pt-rapidity distributions from Kinematics | |
218 | //---------------------------- | |
53399039 | 219 | TH2F *hptyResonanceMC = new TH2F("hptyResonanceMC", " Monte Carlo Resonance",ptbins,ptmin,ptmax,ybins,ymin,ymax); |
220 | hptyResonanceMC->SetLineColor(1); | |
221 | hptyResonanceMC->SetLineStyle(1); | |
222 | hptyResonanceMC->SetLineWidth(2); | |
223 | Ktuple->Project("hptyResonanceMC","y:pt",IdcutResonanceMC); | |
224 | hptyResonanceMC->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
225 | hptyResonanceMC->GetYaxis()->SetTitle("Rapidity"); | |
29e08f93 | 226 | |
227 | cout << " " << endl; | |
228 | cout << "********************" << endl; | |
229 | cout << " Monte Carlo Tracks "<< endl; | |
53399039 | 230 | cout << " " << hptyResonanceMC->GetEntries() << " Resonance in simulation " << endl; |
29e08f93 | 231 | |
232 | ||
233 | //******** Add acceptance cuts - Theta and Pt | |
53399039 | 234 | TH2F *hptyResonanceMCAcc = new TH2F("hptyResonanceMCAcc", " Monte Carlo Resonance in acceptance",ptbins,ptmin,ptmax,ybins,ymin,ymax); |
235 | hptyResonanceMCAcc->SetLineColor(1); | |
236 | hptyResonanceMCAcc->SetLineStyle(1); | |
237 | hptyResonanceMCAcc->SetLineWidth(2); | |
29e08f93 | 238 | |
53399039 | 239 | TString m1MC(IdcutResonanceMC); |
240 | TString m2MC(ResonanceAccCutMC); | |
29e08f93 | 241 | TString m3MC = m1MC + " && " + m2MC ; |
242 | ||
53399039 | 243 | Ktuple->Project("hptyResonanceMCAcc","y:pt",m3MC.Data()); |
244 | hptyResonanceMCAcc->GetYaxis()->SetTitle("Rapidity"); | |
245 | hptyResonanceMCAcc->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
246 | hptyResonanceMCAcc->SetStats(1); | |
247 | hptyResonanceMCAcc->Sumw2(); | |
29e08f93 | 248 | |
53399039 | 249 | TCanvas *c1 = new TCanvas("c1", "Resonance Monte Carlo: Pt vs Y",30,30,700,500); |
29e08f93 | 250 | c1->Range(-1.69394,-0.648855,15.3928,2.77143); |
251 | c1->SetBorderSize(2); | |
252 | c1->SetRightMargin(0.0229885); | |
253 | c1->SetTopMargin(0.0275424); | |
254 | c1->SetFrameFillColor(0); | |
255 | c1->cd(); | |
256 | ||
53399039 | 257 | hptyResonanceMCAcc->SetStats(0); |
258 | hptyResonanceMCAcc->Draw("LEGO2ZFB"); | |
29e08f93 | 259 | //TLatex *tex = new TLatex(2,6,"#mu^{-}"); |
260 | //tex->SetTextSize(0.06); | |
261 | //tex->SetLineWidth(2); | |
262 | //tex->Draw(); | |
263 | ||
53399039 | 264 | sprintf(txt,"Resonance : %d entries",hptyResonanceMCAcc->GetEntries()); |
29e08f93 | 265 | tex = new TLatex(-0.854829,0.794436,txt); |
266 | tex->SetLineWidth(2); | |
267 | tex->Draw(); | |
268 | ||
29e08f93 | 269 | c1->Update(); |
270 | if (SAVE){ | |
53399039 | 271 | c1->SaveAs("ptyResonanceMCAcc.gif"); |
272 | c1->SaveAs("ptyResonanceMCAcc.eps"); | |
29e08f93 | 273 | } |
274 | ||
275 | ||
53399039 | 276 | TH1F *hptResonanceMCAcc = new TH1F("hptResonanceMCAcc", " Monte Carlo Resonance in acceptance",ptbins,ptmin,ptmax); |
277 | hptResonanceMCAcc->SetLineColor(1); | |
278 | hptResonanceMCAcc->SetLineStyle(1); | |
279 | hptResonanceMCAcc->SetLineWidth(2); | |
280 | Ktuple->Project("hptResonanceMCAcc","pt",m3MC.Data()); | |
281 | hptResonanceMCAcc->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
282 | hptResonanceMCAcc->SetStats(1); | |
283 | hptResonanceMCAcc->Sumw2(); | |
29e08f93 | 284 | |
53399039 | 285 | TH1F *hyResonanceMCAcc = new TH1F("hyResonanceMCAcc", " Monte Carlo Resonance in acceptance",ybins,ymin,ymax); |
286 | hyResonanceMCAcc->SetLineColor(1); | |
287 | hyResonanceMCAcc->SetLineStyle(1); | |
288 | hyResonanceMCAcc->SetLineWidth(2); | |
289 | Ktuple->Project("hyResonanceMCAcc","y",m3MC.Data()); | |
290 | hyResonanceMCAcc->GetXaxis()->SetTitle("Rapidity"); | |
291 | hyResonanceMCAcc->SetStats(1); | |
292 | hyResonanceMCAcc->Sumw2(); | |
29e08f93 | 293 | |
294 | ||
295 | // TH2F *hKAccptyCutsMuplus = new TH2F("hKAccptyCutsMuplus", "Monte Carlo #mu^{+} ",ptbins,ptmin,ptmax,ybins,ymin,ymax); | |
296 | // hKAccptyCutsMuplus->SetLineColor(1); | |
297 | // hKAccptyCutsMuplus->SetLineStyle(1); | |
298 | // hKAccptyCutsMuplus->SetLineWidth(2); | |
299 | ||
300 | // TString p1MC(IdcutMuplusMC); | |
301 | // TString p2MC(AccCutMC); | |
302 | // TString p3MC = p1MC + " && " + p2MC ; | |
303 | ||
304 | // Ktuple->Project("hKAccptyCutsMuplus","y:pt",p3MC.Data()); | |
305 | // hKAccptyCutsMuplus->GetYaxis()->SetTitle("Rapidity"); | |
306 | // hKAccptyCutsMuplus->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
307 | // hKAccptyCutsMuplus->SetStats(1); | |
308 | // hKAccptyCutsMuplus->Sumw2(); | |
309 | ||
310 | ||
311 | cout << " " << endl; | |
312 | cout << "********************" << endl; | |
313 | cout << " Monte Carlo Tracks "<< endl; | |
53399039 | 314 | cout << " " << hptyResonanceMCAcc->GetEntries() << " Resonance in acceptance cuts " << endl; |
29e08f93 | 315 | |
316 | ||
317 | /*******************************************************************************************************************/ | |
318 | ||
319 | /*******************************/ | |
320 | /* Reconstructed Tracks Study */ | |
321 | /*******************************/ | |
322 | ||
323 | //---------------------------- | |
324 | // Pt-rapidity distributions from ESD : reconstructed tracks/particle | |
325 | //---------------------------- | |
53399039 | 326 | TH2F *hptyResonanceESD = new TH2F("hptyResonanceESD", " Reconstucted Resonances",ptbins,ptmin,ptmax,ybins,ymin,ymax); |
327 | hptyResonanceESD->SetLineColor(1); | |
328 | hptyResonanceESD->SetLineStyle(1); | |
329 | hptyResonanceESD->SetLineWidth(2); | |
330 | ESDtuple->Project("hptyResonanceESD","y:pt",BckgdCutResonanceESD); | |
331 | hptyResonanceESD->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
332 | hptyResonanceESD->GetYaxis()->SetTitle("Rapidity"); | |
29e08f93 | 333 | |
334 | cout << " " << endl; | |
335 | cout << "********************" << endl; | |
336 | cout << " Reconstructed Tracks" << endl ; | |
53399039 | 337 | cout << " " << hptyResonanceESD->GetEntries() << " Resonance reconstructed " << endl; |
29e08f93 | 338 | |
339 | if (REALISTIC_BACKGROUND){ | |
53399039 | 340 | TH2F *hptyResonanceESDBck = new TH2F("hptyResonanceESDBck", "Background",ptbins,ptmin,ptmax,ybins,ymin,ymax); |
341 | hptyResonanceESDBck->SetLineColor(1); | |
342 | hptyResonanceESDBck->SetLineStyle(1); | |
343 | hptyResonanceESDBck->SetLineWidth(2); | |
344 | ESDtupleBck->Project("hptyResonanceESDBck","y:pt",BckgdCutResonanceESD); | |
345 | hptyResonanceESDBck->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
346 | hptyResonanceESDBck->GetYaxis()->SetTitle("Rapidity"); | |
347 | cout << " with " << hptyResonanceESDBck->GetEntries() << " Resonances from Background (random mixing) " << endl; | |
29e08f93 | 348 | } |
349 | ||
350 | // if something is wrong | |
53399039 | 351 | if ( hptyResonanceESD->GetEntries()==0) {cout << " No entries in hptyResonanceESD " << endl ; break ;} |
29e08f93 | 352 | |
353 | //with Acc cuts - Theta and Pt | |
53399039 | 354 | TH2F *hptyResonanceESDAcc = new TH2F("hptyResonanceESDAcc", "Reconstructed Resonances",ptbins,ptmin,ptmax,ybins,ymin,ymax); |
355 | hptyResonanceESDAcc->SetLineColor(1); | |
356 | hptyResonanceESDAcc->SetLineStyle(1); | |
357 | hptyResonanceESDAcc->SetLineWidth(2); | |
29e08f93 | 358 | |
53399039 | 359 | TString m1Rec(BckgdCutResonanceESD); |
360 | TString m2Rec(ResonanceAccCutESD); | |
29e08f93 | 361 | TString m3Rec = m1Rec + " && " + m2Rec ; |
362 | ||
53399039 | 363 | ESDtuple->Project("hptyResonanceESDAcc","y:pt",m3Rec.Data()); |
364 | hptyResonanceESDAcc->GetYaxis()->SetTitle("Rapidity"); | |
365 | hptyResonanceESDAcc->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
366 | hptyResonanceESDAcc->SetStats(1); | |
367 | hptyResonanceESDAcc->Sumw2(); | |
29e08f93 | 368 | |
369 | cout << " " << endl; | |
370 | cout << "********************" << endl; | |
371 | cout << " Reconstructed Tracks" << endl ; | |
53399039 | 372 | cout << " " << hptyResonanceESDAcc->GetEntries() << " Resonance in acceptance cuts " << endl; |
29e08f93 | 373 | |
374 | if (REALISTIC_BACKGROUND){ | |
375 | //with Acc cuts - Theta and Pt | |
53399039 | 376 | TH2F *hptyResonanceESDBckAcc = new TH2F("hptyResonanceESDBckAcc", "Reconstructed Resonances",ptbins,ptmin,ptmax,ybins,ymin,ymax); |
377 | hptyResonanceESDBckAcc->SetLineColor(1); | |
378 | hptyResonanceESDBckAcc->SetLineStyle(1); | |
379 | hptyResonanceESDBckAcc->SetLineWidth(2); | |
29e08f93 | 380 | |
53399039 | 381 | TString m1Rec(BckgdCutResonanceESD); |
382 | TString m2Rec(ResonanceAccCutESD); | |
29e08f93 | 383 | TString m3Rec = m1Rec + " && " + m2Rec ; |
384 | ||
53399039 | 385 | ESDtupleBck->Project("hptyResonanceESDBckAcc","y:pt",m3Rec.Data()); |
386 | hptyResonanceESDBckAcc->GetYaxis()->SetTitle("Rapidity"); | |
387 | hptyResonanceESDBckAcc->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
388 | hptyResonanceESDBckAcc->SetStats(1); | |
389 | hptyResonanceESDBckAcc->Sumw2(); | |
390 | cout << " with " << hptyResonanceESDBckAcc->GetEntries() << " Resonances from Background (random mixing) " << endl; | |
29e08f93 | 391 | } |
392 | ||
393 | ||
53399039 | 394 | TCanvas *c100 = new TCanvas("c100", "Resonance Reconstructed in Acceptance: Pt vs Y",210,30,700,500); |
29e08f93 | 395 | c100->Range(-1.69394,-0.648855,15.3928,2.77143); |
396 | c100->SetBorderSize(2); | |
397 | c100->SetRightMargin(0.0229885); | |
398 | c100->SetTopMargin(0.0275424); | |
399 | c100->SetFrameFillColor(0); | |
400 | ||
401 | c100->cd(); | |
53399039 | 402 | hptyResonanceESDAcc->SetStats(0); |
403 | hptyResonanceESDAcc->Draw("LEGO2ZFB"); | |
404 | sprintf(txt,"Resonance : %d entries",hptyResonanceESDAcc->GetEntries()); | |
29e08f93 | 405 | tex = new TLatex(-0.854829,0.794436,txt); |
406 | tex->SetLineWidth(2); | |
407 | tex->Draw(); | |
408 | ||
409 | c100->Update(); | |
410 | if (SAVE){ | |
53399039 | 411 | c100->SaveAs("ptyResonanceESDAcc.gif"); |
412 | c100->SaveAs("ptyResonanceESDAcc.eps"); | |
29e08f93 | 413 | } |
414 | ||
29e08f93 | 415 | if (REALISTIC_BACKGROUND){ |
53399039 | 416 | TCanvas *c110 = new TCanvas("c110", "Resonance Background Reconstructed in Acceptance: Pt vs Y",215,35,700,500); |
29e08f93 | 417 | c110->Range(-1.69394,-0.648855,15.3928,2.77143); |
418 | c110->SetBorderSize(2); | |
419 | c110->SetRightMargin(0.0229885); | |
420 | c110->SetTopMargin(0.0275424); | |
421 | c110->SetFrameFillColor(0); | |
422 | ||
423 | c110->cd(); | |
53399039 | 424 | hptyResonanceESDBckAcc->SetStats(0); |
425 | hptyResonanceESDBckAcc->Draw("LEGO2ZFB"); | |
426 | sprintf(txt,"Resonance backround : %d entries",hptyResonanceESDBckAcc->GetEntries()); | |
29e08f93 | 427 | tex = new TLatex(-0.854829,0.794436,txt); |
428 | tex->SetLineWidth(2); | |
429 | tex->Draw(); | |
430 | ||
431 | ||
432 | c110->Update(); | |
433 | if (SAVE){ | |
53399039 | 434 | c110->SaveAs("ptyResonanceESDBckAcc.gif"); |
435 | c110->SaveAs("ptyResonanceESDBckAcc.eps"); | |
29e08f93 | 436 | } |
437 | } | |
438 | ||
53399039 | 439 | TH1F *hptResonanceESDAcc = new TH1F("hptResonanceESDAcc", " Monte Carlo Resonance in acceptance",ptbins,ptmin,ptmax); |
440 | hptResonanceESDAcc->SetLineColor(1); | |
441 | hptResonanceESDAcc->SetLineStyle(1); | |
442 | hptResonanceESDAcc->SetLineWidth(2); | |
443 | ESDtuple->Project("hptResonanceESDAcc","pt",m3Rec.Data()); | |
444 | hptResonanceESDAcc->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
445 | hptResonanceESDAcc->SetStats(1); | |
446 | hptResonanceESDAcc->Sumw2(); | |
447 | ||
448 | TH1F *hyResonanceESDAcc = new TH1F("hyResonanceESDAcc", " Monte Carlo Resonance in acceptance",ybins,ymin,ymax); | |
449 | hyResonanceESDAcc->SetLineColor(1); | |
450 | hyResonanceESDAcc->SetLineStyle(1); | |
451 | hyResonanceESDAcc->SetLineWidth(2); | |
452 | ESDtuple->Project("hyResonanceESDAcc","y",m3Rec.Data()); | |
453 | hyResonanceESDAcc->GetXaxis()->SetTitle("Rapidity"); | |
454 | hyResonanceESDAcc->SetStats(1); | |
455 | hyResonanceESDAcc->Sumw2(); | |
29e08f93 | 456 | |
457 | ||
458 | /*******************************************************************************************************************/ | |
459 | ||
460 | /*******************************/ | |
461 | /* Efficiencies calculations */ | |
462 | /*******************************/ | |
463 | ||
464 | ||
465 | cout << " " << endl; | |
466 | cout << "***********************" << endl; | |
467 | cout << " Integrated efficiency" << endl ; | |
468 | ||
469 | if (REALISTIC_BACKGROUND) | |
53399039 | 470 | cout << " ResonanceESDAcc/ResonanceMCAcc = " << (hptyResonanceESDAcc->GetEntries()-hptyResonanceESDBckAcc->GetEntries())/hptyResonanceMCAcc->GetEntries() << endl; |
29e08f93 | 471 | else |
53399039 | 472 | cout << " ResonanceESDAcc/ResonanceMCAcc = " << hptyResonanceESDAcc->GetEntries()/hptyResonanceMCAcc->GetEntries() << endl; |
29e08f93 | 473 | |
53399039 | 474 | TH2F *hEffptyResonance = new TH2F("hEffptyResonance", " Resonance Efficiency",ptbins,ptmin,ptmax,ybins,ymin,ymax); |
29e08f93 | 475 | |
476 | if (REALISTIC_BACKGROUND){ | |
53399039 | 477 | TH2F *hptyResonanceESDAccBckSubstracted = new TH2F("hptyResonanceESDAccBckSubstracted","hptyResonanceESDAccBckSubstracted",ptbins,ptmin,ptmax,ybins,ymin,ymax); |
478 | hptyResonanceESDAccBckSubstracted->Add(hptyResonanceESDAcc,hptyResonanceESDBckAcc,1,-1); | |
479 | hEffptyResonance->Divide(hptyResonanceESDAccBckSubstracted,hptyResonanceMCAcc,1,1); | |
29e08f93 | 480 | } |
481 | else | |
53399039 | 482 | hEffptyResonance->Divide(hptyResonanceESDAcc,hptyResonanceMCAcc,1,1); |
29e08f93 | 483 | |
53399039 | 484 | TCanvas *c1000 = new TCanvas("c1000", "Resonance efficiency : Pt vs Y",390,30,700,500); |
29e08f93 | 485 | c1000->Range(-1.69394,-0.648855,15.3928,2.77143); |
486 | c1000->SetBorderSize(2); | |
487 | c1000->SetRightMargin(0.0229885); | |
488 | c1000->SetTopMargin(0.0275424); | |
489 | c1000->SetFrameFillColor(0); | |
490 | ||
491 | c1000->cd(); | |
53399039 | 492 | hEffptyResonance->SetStats(0); |
493 | hEffptyResonance->Draw("LEGO2fz"); | |
29e08f93 | 494 | |
495 | ||
496 | ||
53399039 | 497 | TH1F *hEffptResonance = new TH1F("hEffptResonance", "Resonance Efficiency vs pt",ptbins,ptmin,ptmax); |
498 | hEffptResonance->Divide(hptResonanceESDAcc,hptResonanceMCAcc,1,1); | |
499 | hEffptResonance->SetLineWidth(2); | |
500 | hEffptResonance->SetMinimum(0); | |
501 | hEffptResonance->SetStats(1); | |
29e08f93 | 502 | |
503 | ||
53399039 | 504 | TCanvas *c1100 = new TCanvas("c1100", "Resonance efficiency : Pt",410,50,700,500); |
29e08f93 | 505 | c1100->Range(-1.69394,-0.648855,15.3928,2.77143); |
506 | c1100->SetBorderSize(2); | |
507 | c1100->SetRightMargin(0.0229885); | |
508 | c1100->SetTopMargin(0.0275424); | |
509 | c1100->SetFrameFillColor(0); | |
510 | ||
511 | c1100->cd(); | |
53399039 | 512 | hEffptResonance->SetStats(0); |
513 | hEffptResonance->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
514 | hEffptResonance->GetYaxis()->SetTitle("Efficiency "); | |
515 | hEffptResonance->Draw("E"); | |
29e08f93 | 516 | |
517 | ||
518 | ||
53399039 | 519 | TH1F *hEffyResonance = new TH1F("hEffyResonance", "Resonance Efficiency vs y",ybins,ymin,ymax); |
520 | hEffyResonance->Divide(hyResonanceESDAcc,hyResonanceMCAcc,1,1); | |
521 | hEffyResonance->SetLineWidth(2); | |
522 | hEffyResonance->SetStats(1); | |
29e08f93 | 523 | |
53399039 | 524 | TCanvas *c1200 = new TCanvas("c1200", "Resonance efficiency : Y",430,70,700,500); |
29e08f93 | 525 | c1200->Range(-1.69394,-0.648855,15.3928,2.77143); |
526 | c1200->SetBorderSize(2); | |
527 | c1200->SetRightMargin(0.0229885); | |
528 | c1200->SetTopMargin(0.0275424); | |
529 | c1200->SetFrameFillColor(0); | |
530 | ||
531 | c1200->cd(); | |
53399039 | 532 | hEffyResonance->SetStats(0); |
533 | hEffyResonance->GetXaxis()->SetTitle("Rapidity"); | |
534 | hEffyResonance->GetYaxis()->SetTitle("Efficiency "); | |
535 | hEffyResonance->Draw("E"); | |
29e08f93 | 536 | |
537 | c1000->Update(); | |
538 | c1100->Update(); | |
539 | c1200->Update(); | |
540 | if (SAVE){ | |
53399039 | 541 | c1000->SaveAs("EffptyResonance.gif"); |
542 | c1000->SaveAs("EffptyResonance.eps"); | |
543 | c1100->SaveAs("EffptResonance.gif"); | |
544 | c1100->SaveAs("EffptResonance.eps"); | |
545 | c1200->SaveAs("EffyResonance.gif"); | |
546 | c1200->SaveAs("EffyResonance.eps"); | |
29e08f93 | 547 | } |
548 | ||
549 | /*******************************************************************************************************************/ | |
550 | ||
551 | /*******************************/ | |
552 | /* Trigger matching */ | |
553 | /*******************************/ | |
554 | ||
555 | ||
556 | Float_t triggerChi2Min = 0.; | |
557 | Float_t triggerChi2Max = 7.5; | |
558 | ||
53399039 | 559 | //TString m1Rec(BckgdCutResonanceESD); |
560 | //TString m2Rec(ResonanceAccCutESD); | |
29e08f93 | 561 | TString m1TG = m1Rec + " && " + m2Rec ; |
562 | //cout << m1TG.Data() << endl; | |
563 | ||
564 | TH1F *hInvMassNoTriggerCut= new TH1F("hInvMassNoTriggerCut","hInvMassNoTriggerCut",invMassBins,invMassMin,invMassMax); | |
565 | hInvMassNoTriggerCut->GetXaxis()->SetTitle("Mass [GeV/c^{2}]"); | |
566 | hInvMassNoTriggerCut->GetYaxis()->SetTitle("Counts"); | |
567 | ESDtuple->Project("hInvMassNoTriggerCut","minv",m1TG.Data()); | |
568 | ||
569 | //Add trigger UnlikePairAllPt | |
570 | TString m2TG = m1TG + " && (tw & 0x800) == 2048" ; | |
53399039 | 571 | TH1F *hInvMassResonanceTrigger= new TH1F("hInvMassResonanceTrigger","hInvMassResonanceTrigger",invMassBins,invMassMin,invMassMax); |
572 | hInvMassResonanceTrigger->GetXaxis()->SetTitle("Mass [GeV/c^{2}]"); | |
573 | hInvMassResonanceTrigger->GetYaxis()->SetTitle("Counts"); | |
574 | ESDtuple->Project("hInvMassResonanceTrigger","minv",m2TG.Data()); | |
29e08f93 | 575 | |
576 | ||
577 | cout << " " << endl; | |
578 | cout << "********************" << endl; | |
579 | cout << " TriggerMatching" << endl ; | |
53399039 | 580 | cout << " " << hInvMassResonanceTrigger->GetEntries() << " Resonance with trigger UnlikePairAllPt " << endl; |
29e08f93 | 581 | |
582 | TString m2TGNo = m1TG + " && (tw & 0x800) != 2048" ; | |
53399039 | 583 | TH1F *hResonanceTriggerNo= new TH1F("hResonanceTriggerNo","hResonanceTriggerNo",32769,0,32768); |
584 | hResonanceTriggerNo->GetXaxis()->SetTitle("Mass [GeV/c^{2}]"); | |
585 | hResonanceTriggerNo->GetYaxis()->SetTitle("Counts"); | |
586 | ESDtuple->Project("hResonanceTriggerNo","tw",m2TGNo.Data()); | |
29e08f93 | 587 | |
588 | ||
589 | ||
590 | ||
591 | //Add matching rec/trig for 2 tracks | |
592 | TString m3TG = m2TG + " && trig1 > 0 && trig2 > 0" ; | |
53399039 | 593 | TH1F *hInvMassResonanceTriggerTwoMatch = new TH1F("hInvMassResonanceTriggerTwoMatch","hInvMassResonanceTrigger with 2 matching tracks",invMassBins,invMassMin,invMassMax); |
594 | hInvMassResonanceTriggerTwoMatch->GetXaxis()->SetTitle("Mass [GeV/c^{2}]"); | |
595 | hInvMassResonanceTriggerTwoMatch->GetYaxis()->SetTitle("Counts"); | |
596 | ESDtuple->Project("hInvMassResonanceTriggerTwoMatch","minv",m3TG.Data()); | |
29e08f93 | 597 | |
598 | ||
599 | //Add matching rec/trig for 1 tracks | |
600 | TString m4TG = m2TG + " && (trig1 > 0 || trig2 > 0)" ; | |
53399039 | 601 | TH1F *hInvMassResonanceTriggerOneMatch= new TH1F("hInvMassResonanceTriggerOneMatch","hInvMassResonanceTrigger with 1 matching tracks",invMassBins,invMassMin,invMassMax); |
602 | hInvMassResonanceTriggerOneMatch->GetXaxis()->SetTitle("Mass [GeV/c^{2}]"); | |
603 | hInvMassResonanceTriggerOneMatch->GetYaxis()->SetTitle("Counts"); | |
604 | ESDtuple->Project("hInvMassResonanceTriggerOneMatch","minv",m4TG.Data()); | |
29e08f93 | 605 | |
606 | TString m5TG = m2TG + " && (trig1 == 0 && trig2 == 0)" ; | |
53399039 | 607 | TH1F *hInvMassResonanceTriggerNoMatch= new TH1F("hInvMassResonanceTriggerNoMatch","hInvMassResonanceTrigger with no matching tracks",invMassBins,invMassMin,invMassMax); |
608 | hInvMassResonanceTriggerNoMatch->GetXaxis()->SetTitle("Mass [GeV/c^{2}]"); | |
609 | hInvMassResonanceTriggerNoMatch->GetYaxis()->SetTitle("Counts"); | |
610 | ESDtuple->Project("hInvMassResonanceTriggerNoMatch","minv",m5TG.Data()); | |
29e08f93 | 611 | |
53399039 | 612 | TCanvas *c2 = new TCanvas("c2", "Resonance Trigger efficiency",30,90,700,500); |
29e08f93 | 613 | c2->Range(-1.69394,-0.648855,15.3928,2.77143); |
614 | c2->SetBorderSize(2); | |
615 | c2->SetRightMargin(0.0229885); | |
616 | c2->SetTopMargin(0.0275424); | |
617 | c2->SetFrameFillColor(0); | |
618 | c2->SetLogy(1); | |
619 | ||
620 | c2->cd(); | |
621 | hInvMassNoTriggerCut->SetStats(0); | |
622 | hInvMassNoTriggerCut->SetLineWidth(2); | |
623 | hInvMassNoTriggerCut->SetLineStyle(3); | |
624 | hInvMassNoTriggerCut->SetLineColor(1); | |
625 | hInvMassNoTriggerCut->GetYaxis()->SetTitleOffset(1.2); | |
626 | hInvMassNoTriggerCut->Draw(); | |
53399039 | 627 | hInvMassResonanceTrigger->SetLineWidth(2); |
628 | hInvMassResonanceTrigger->SetLineStyle(0); | |
629 | hInvMassResonanceTrigger->SetLineColor(2); | |
630 | hInvMassResonanceTrigger->Draw("same"); | |
631 | hInvMassResonanceTriggerTwoMatch->SetLineWidth(2); | |
632 | hInvMassResonanceTriggerTwoMatch->SetLineStyle(1); | |
633 | hInvMassResonanceTriggerTwoMatch->SetLineColor(4); | |
634 | hInvMassResonanceTriggerTwoMatch->Draw("same"); | |
635 | hInvMassResonanceTriggerOneMatch->SetLineWidth(2); | |
636 | hInvMassResonanceTriggerOneMatch->SetLineStyle(2); | |
637 | hInvMassResonanceTriggerOneMatch->SetLineColor(51); | |
638 | hInvMassResonanceTriggerOneMatch->Draw("same"); | |
639 | hInvMassResonanceTriggerNoMatch->SetLineWidth(2); | |
640 | hInvMassResonanceTriggerNoMatch->SetLineStyle(2); | |
641 | hInvMassResonanceTriggerNoMatch->SetLineColor(46); | |
642 | hInvMassResonanceTriggerNoMatch->Draw("same"); | |
29e08f93 | 643 | |
644 | TLegend *leg = new TLegend(0.12,0.6,0.50,0.89); | |
53399039 | 645 | leg->SetHeader("Reconstructed Resonance Invariant Mass"); |
29e08f93 | 646 | |
647 | leg->AddEntry(hInvMassNoTriggerCut,Form("All (%.0f cnts)",hInvMassNoTriggerCut->GetEntries()),"l"); | |
53399039 | 648 | leg->AddEntry(hInvMassResonanceTrigger,Form("UnlikePairAllPt Trigger (%.0f cnts)",hInvMassResonanceTrigger->GetEntries()),"l"); |
649 | leg->AddEntry(hInvMassResonanceTriggerTwoMatch,Form("UPAllPt Trig. and 2 tracks matches (%.0f cnts)",hInvMassResonanceTriggerTwoMatch->GetEntries()),"l"); | |
650 | leg->AddEntry(hInvMassResonanceTriggerOneMatch,Form("UPAllPt Trig. and 1 track match (%.0f cnts)",hInvMassResonanceTriggerOneMatch->GetEntries()),"l"); | |
651 | leg->AddEntry(hInvMassResonanceTriggerNoMatch,Form("UPAllPt Trig. and no matching track (%.0f cnts)",hInvMassResonanceTriggerNoMatch->GetEntries()),"l"); | |
29e08f93 | 652 | leg->Draw(); |
653 | ||
654 | c2->Update(); | |
655 | if(SAVE){ | |
656 | c2->SaveAs("TriggerMatching.gif"); | |
657 | c2->SaveAs("TriggerMatching.eps"); | |
658 | } | |
659 | ||
660 | /*******************************************************************************************************************/ | |
661 | ||
662 | /*******************************/ | |
663 | /* Mass resolution */ | |
664 | /*******************************/ | |
665 | ||
666 | const Int_t nofMassHistograms = ptbins ; | |
53399039 | 667 | // cs invMassMin = 7.0 ; |
668 | // cs invMassMax = 11.0 ; | |
669 | // cs invMassBins = 100 ; | |
29e08f93 | 670 | |
53399039 | 671 | //Float_t ptmin = 0.0; |
29e08f93 | 672 | //Float_t ptmax = 16.0; Int_t ptbins = 8; |
673 | Float_t ptbinssize = ptmax/ptbins; | |
53399039 | 674 | |
29e08f93 | 675 | |
676 | Float_t ptbinslimits[nofMassHistograms+1]; | |
677 | Float_t ptbinscenters[nofMassHistograms]; | |
678 | ||
679 | for(Int_t as = 0 ; as < nofMassHistograms+1 ; as++){ | |
680 | ptbinslimits[as] = as*ptbinssize ; | |
681 | //cout << "ptbinslimits["<<as<<"] = " << ptbinslimits[as] << endl ; | |
682 | } | |
683 | for(Int_t as = 0 ; as < nofMassHistograms ; as++) { | |
684 | ptbinscenters[as] = ptbinslimits[as] + (ptbinslimits[as+1] - ptbinslimits[as])/2 ; | |
685 | //cout << "ptbinscenters[as] = " << ptbinscenters[as] << endl ; | |
686 | } | |
687 | ||
688 | ||
689 | TH1F **hInvMassInPtBins = new TH1F* [nofMassHistograms] ; | |
690 | TString histExt; | |
691 | Char_t theExt[20]; | |
692 | TString histCut; | |
693 | Char_t theCut[100]; | |
694 | ||
695 | ||
696 | // Load the fit functions | |
697 | //gROOT->ProcessLine(".L /home/suire/AliRoot/macros/FIT/FitFunction.C"); | |
698 | ||
699 | TF1 *fitFunc ; | |
700 | ||
53399039 | 701 | if (fitfunc==0 || fitfunc==1){ |
702 | fitFunc = new TF1("gaussian","gaus(0)",FitLow,FitHigh); | |
703 | if (!fitFunc) | |
704 | fitFunc = new TF1("gaussian","[0]*TMath::Exp(-0.5*(x-[1])*(x-[1])/[2]/[2])",FitLow,FitHigh); | |
705 | fitFunc->SetParNames("Constant","Mean value","Width"); | |
29e08f93 | 706 | fitFunc->SetLineColor(2); |
707 | fitFunc->SetLineWidth(2); | |
708 | fitFunc->SetLineStyle(2); | |
709 | ||
53399039 | 710 | Float_t *tParams = new Float_t[3] ; |
711 | if (ResType==553){ | |
712 | tParams[0] = 500; | |
713 | tParams[1] = 9.47; | |
714 | tParams[2] = 0.1; | |
715 | } | |
716 | if (ResType==443){ | |
717 | tParams[0] = 500; | |
718 | tParams[1] = 3.1; | |
719 | tParams[2] = 0.1; | |
720 | } | |
29e08f93 | 721 | |
53399039 | 722 | fitFunc->SetParameters(tParams[0],tParams[1],tParams[2]); |
29e08f93 | 723 | |
53399039 | 724 | TString FitFuncName = "gaussian"; |
29e08f93 | 725 | } |
726 | ||
727 | if (fitfunc==2){ | |
728 | fitFunc = new TF1("fitlan_a",fitlan_a,FitLow,FitHigh,3); | |
729 | fitFunc->SetParNames("Constant","Mean value","Width"); | |
730 | fitFunc->SetLineColor(2); | |
731 | fitFunc->SetLineWidth(2); | |
732 | fitFunc->SetLineStyle(2); | |
733 | ||
734 | Float_t *tParams = new Float_t[3] ; | |
735 | tParams[0] = 1000; | |
736 | tParams[1] = 9.47; | |
737 | tParams[2] = 0.05; | |
738 | fitFunc->SetParameters(tParams[0],tParams[1],tParams[2]); | |
739 | ||
740 | TString FitFuncName = "fitlan_a"; | |
741 | } | |
742 | ||
53399039 | 743 | if (fitfunc==3){ |
744 | fitFunc = new TF1("fitlangaussUpsilon",fitlangaussUpsilon,FitLow,FitHigh,6); | |
745 | fitFunc->SetParNames("Constant","Mean value","Width","SigmaGauss","LowFitVal","HighFitVal"); | |
29e08f93 | 746 | fitFunc->SetLineColor(2); |
747 | fitFunc->SetLineWidth(2); | |
748 | fitFunc->SetLineStyle(2); | |
749 | ||
53399039 | 750 | Float_t *tParams = new Float_t[6] ; |
751 | tParams[0] = 5000; | |
752 | tParams[1] = 9.47; | |
753 | tParams[2] = 0.05; //0.5 | |
754 | tParams[3] = 0.05; // 1. | |
29e08f93 | 755 | |
53399039 | 756 | tParams[4] = FitLow; |
757 | tParams[5] = FitHigh; | |
758 | ||
759 | fitFunc->SetParameters(tParams[0],tParams[1],tParams[2],tParams[3],tParams[4],tParams[5]); | |
760 | fitFunc->FixParameter(4,FitLow); | |
761 | fitFunc->FixParameter(5,FitHigh); | |
762 | ||
763 | fitFunc->SetParLimits(1, 9.1 ,9.7); | |
764 | fitFunc->SetParLimits(2, 0.001 ,0.1); | |
765 | fitFunc->SetParLimits(3, 0.001 ,0.5); | |
766 | ||
767 | // for special case one may use | |
768 | //fitFunc->SetParameter(1,9.35); | |
769 | ||
770 | TString FitFuncName = "fitlangaussUpsilon"; | |
29e08f93 | 771 | } |
772 | ||
773 | if (fitfunc==4){ | |
774 | fitFunc = new TF1("fitlanUpsilon",fitlanUpsilon,FitLow,FitHigh,3); | |
775 | fitFunc->SetParNames("Constant","Mean value","Width"); | |
776 | fitFunc->SetLineColor(2); | |
777 | fitFunc->SetLineWidth(2); | |
778 | fitFunc->SetLineStyle(2); | |
779 | ||
780 | Float_t *tParams = new Float_t[3] ; | |
781 | tParams[0] = 1000; | |
782 | tParams[1] = 9.47; | |
783 | tParams[2] = 0.05; | |
784 | fitFunc->SetParameters(tParams[0],tParams[1],tParams[2]); | |
785 | ||
786 | TString FitFuncName = "fitlanUpsilon"; | |
787 | } | |
788 | ||
789 | ||
790 | Float_t fitResultsWidth[nofMassHistograms]; | |
791 | Float_t fitResultsWidthErr[nofMassHistograms]; | |
792 | ||
793 | Float_t fitResultsMean[nofMassHistograms]; | |
794 | Float_t fitResultsMeanErr[nofMassHistograms]; | |
795 | ||
796 | Float_t errx[nofMassHistograms]; | |
797 | ||
798 | ||
799 | ||
53399039 | 800 | TCanvas *call = new TCanvas("call", "Resonance : invariant mass spectra",30,330,700,500); |
29e08f93 | 801 | call->Range(-1.69394,-0.648855,15.3928,2.77143); |
802 | call->SetBorderSize(2); | |
803 | //call->SetRightMargin(0.2229885); | |
804 | call->SetRightMargin(0.0229885); | |
805 | call->SetTopMargin(0.0275424); | |
806 | call->SetFrameFillColor(0); | |
807 | call->cd(); | |
53399039 | 808 | |
809 | TH1F* hInvMass = new TH1F("hInvMass","Inv. Mass Pt,Y integrated",30,0,3+RESONANCE_MASS); | |
29e08f93 | 810 | ESDtuple->Project("hInvMass","minv"); |
811 | hInvMass->SetLineWidth(2); | |
812 | hInvMass->Draw("HE"); | |
53399039 | 813 | |
29e08f93 | 814 | if(REALISTIC_BACKGROUND){ |
53399039 | 815 | TH1F* hInvMassBck = new TH1F("hInvMassBck","Background Pt,Y integrated",30,0,3+RESONANCE_MASS); |
29e08f93 | 816 | ESDtupleBck->Project("hInvMassBck","minv"); |
817 | hInvMassBck->SetLineWidth(2); | |
818 | hInvMassBck->SetLineStyle(2); | |
819 | hInvMassBck->SetLineColor(2); | |
820 | hInvMassBck->Draw("same"); | |
821 | } | |
822 | call->Modified(); | |
823 | call->Update(); | |
824 | ||
825 | ||
53399039 | 826 | TCanvas *cfit = new TCanvas("cfit", "Resonance : invariant mass fit",30,330,700,500); |
29e08f93 | 827 | cfit->Range(-1.69394,-0.648855,15.3928,2.77143); |
828 | cfit->SetBorderSize(2); | |
829 | //cfit->SetRightMargin(0.2229885); | |
830 | cfit->SetRightMargin(0.0229885); | |
831 | cfit->SetTopMargin(0.0275424); | |
832 | cfit->SetFrameFillColor(0); | |
833 | ||
834 | //TH1F *hInvMassAllClone = (TH1F*) gROOT->FindObject("hInvMassAll"); | |
835 | ||
836 | TH1F* hInvMassAll = new TH1F("hInvMassAll","Inv. Mass Pt,Y integrated",invMassBins,invMassMin,invMassMax); | |
837 | sprintf(theCut,"pt> %.2f && pt<= %.2f",ptbinslimits[0],ptbinslimits[nofMassHistograms]); | |
838 | //cout << "theCut" << theCut << endl ; | |
839 | ESDtuple->Project("hInvMassAll","minv",theCut); | |
840 | hInvMassAll->Draw(); | |
841 | ||
842 | cfit->Update(); | |
843 | if (SAVE){ | |
53399039 | 844 | cfit->SaveAs("ResonanceMass.gif"); |
845 | cfit->SaveAs("ResonanceMass.eps"); | |
29e08f93 | 846 | } |
847 | ||
848 | if(REALISTIC_BACKGROUND){ | |
849 | TH1F* hInvMassAllBck = new TH1F("hInvMassAllBck","Background Pt,Y integrated",invMassBins,invMassMin,invMassMax); | |
850 | ESDtupleBck->Project("hInvMassAllBck","minv",theCut); | |
851 | hInvMassAllBck->SetLineWidth(2); | |
852 | hInvMassAllBck->SetLineStyle(2); | |
853 | hInvMassAllBck->SetLineColor(2); | |
854 | hInvMassAllBck->Draw("same"); | |
855 | } | |
856 | cfit->Modified(); | |
857 | cfit->Update(); | |
858 | ||
29e08f93 | 859 | //Fit also the pt-integrated mass histogram |
860 | cfit->cd(); | |
861 | hInvMassAll->Fit(FitFuncName.Data(),"rv"); | |
862 | if (FitFuncName == "gaussian") | |
863 | fitFunc->SetRange(fitFunc->GetParameter(1)-gaussianFitWidth,fitFunc->GetParameter(1)+gaussianFitWidth); | |
864 | hInvMassAll->Fit(FitFuncName.Data(),"rv"); | |
865 | cfit->Modified(); | |
866 | cfit->Update(); | |
867 | ||
868 | cout << " " << endl; | |
869 | cout << "********************************************" << endl; | |
870 | cout << " Fit ("<<FitFuncName.Data()<<" results of the pt-integrated mass peak " << endl ; | |
871 | cout << " Mean value : " << fitFunc->GetParameter(1) << " +/- " << fitFunc->GetParError(1) << endl ; | |
872 | cout << " Width value : " << fitFunc->GetParameter(2) << " +/- " << fitFunc->GetParError(2) << endl ; | |
873 | cout << " ChiSquare of the fit : " << fitFunc->GetChisquare()<< " / " << fitFunc->GetNDF() << endl ; | |
874 | cout << "********************************************" << endl; | |
875 | cout << " " << endl; | |
29e08f93 | 876 | |
53399039 | 877 | TCanvas *cptfits = new TCanvas("cptfits", "Resonance : invariant mass fit",30,330,700,500); |
29e08f93 | 878 | cptfits->Range(-1.69394,-0.648855,15.3928,2.77143); |
879 | cptfits->SetBorderSize(2); | |
880 | //cptfits->SetRightMargin(0.2229885); | |
881 | cptfits->SetRightMargin(0.0229885); | |
882 | cptfits->SetTopMargin(0.0275424); | |
883 | cptfits->SetFrameFillColor(0); | |
884 | ||
885 | Float_t chi2sum = 0.; | |
886 | ||
887 | for (Int_t qw = 0 ; qw < nofMassHistograms ; qw++){ | |
888 | sprintf(theExt,"_%d_%d",ptbinslimits[qw],ptbinslimits[qw+1]); | |
889 | histExt= theExt; | |
890 | hInvMassInPtBins[qw] = new TH1F("hInvMassInPtBins"+histExt,"hInvMassInPtBins"+histExt,invMassBins,invMassMin,invMassMax); | |
891 | hInvMassInPtBins[qw]->GetXaxis()->SetTitle("Mass [GeV/c^{2}]"); | |
892 | hInvMassInPtBins[qw]->GetYaxis()->SetTitle("Counts"); | |
893 | ||
894 | sprintf(theCut,"pt> %.2f && pt<= %.2f",ptbinslimits[qw],ptbinslimits[qw+1]); | |
895 | histCut = theCut; | |
896 | ESDtuple->Project("hInvMassInPtBins"+histExt,"minv",histCut.Data()); | |
897 | ||
898 | if (hInvMassInPtBins[qw]->GetEntries() > 80){ | |
899 | cptfits->cd(); | |
900 | ||
901 | if (FitFuncName == "gaussian"){ | |
902 | //prefit with a gaussian to get the maximum set a correct range | |
903 | hInvMassInPtBins[qw]->Fit(FitFuncName.Data(),"rv"); | |
904 | fitFunc->SetRange(fitFunc->GetParameter(1)-gaussianFitWidth,fitFunc->GetParameter(1)+gaussianFitWidth); | |
905 | } | |
906 | ||
907 | hInvMassInPtBins[qw]->SetStats(1); | |
908 | hInvMassInPtBins[qw]->Draw("HE"); | |
909 | hInvMassInPtBins[qw]->Fit(FitFuncName.Data(),"rv"); | |
910 | ||
911 | ||
912 | cout << "==============" << endl; | |
913 | cout << "ChiSquare of the fit : " << fitFunc->GetChisquare()<< " / " << fitFunc->GetNDF() | |
914 | << " = " << fitFunc->GetChisquare()/fitFunc->GetNDF() <<endl; | |
915 | cout << "==============" << endl; | |
916 | chi2sum += fitFunc->GetChisquare()/fitFunc->GetNDF(); | |
917 | ||
918 | cptfits->Modified(); | |
919 | cptfits->Update(); | |
920 | ||
921 | fitResultsWidth[qw] = TMath::Abs(fitFunc->GetParameter(2)); | |
922 | fitResultsWidthErr[qw] = fitFunc->GetParError(2); | |
923 | ||
924 | fitResultsMean[qw] = TMath::Abs(fitFunc->GetParameter(1)); | |
925 | fitResultsMeanErr[qw] = fitFunc->GetParError(1); | |
926 | ||
927 | if(FitFuncName == "fitlangaussUpsilon"){ | |
928 | // width = gauss_width + width landau | |
929 | fitResultsWidth[qw] = TMath::Abs(fitFunc->GetParameter(2)) + TMath::Abs(fitFunc->GetParameter(3)); | |
930 | fitResultsWidthErr[qw] = TMath::Abs(fitFunc->GetParError(2)) + TMath::Abs(fitFunc->GetParError(3)); | |
931 | ||
932 | // problem with the mean of the fit (parameter1) which doesn't give the peak position | |
933 | fitResultsMean[qw] = TMath::Abs(fitFunc->GetMaximumX()); | |
934 | fitResultsMeanErr[qw] = fitFunc->GetParError(1); | |
935 | //this could work but it's not very nice.... | |
936 | } | |
937 | ||
938 | } | |
939 | else { | |
940 | fitResultsWidth[qw] = 0.; | |
941 | fitResultsWidthErr[qw] = 0.; | |
942 | fitResultsMean[qw] = 0.; | |
943 | fitResultsMeanErr[qw] = 0.; | |
944 | } | |
945 | ||
946 | errx[qw] = 0.; | |
947 | ||
948 | } | |
949 | ||
950 | ||
951 | //From fit results extract a Mean value for the Mass and width | |
952 | cout << " " << endl; | |
953 | cout << "********************************************" << endl; | |
954 | cout << " Fit Function : " << FitFuncName << " " ; | |
955 | if (FitFuncName == "gaussian" ) cout << "+/- " << gaussianFitWidth << " MeV/c2" ; | |
956 | cout << endl ; | |
957 | cout << " Fit results : mean values " << endl ; | |
958 | Float_t meanMass = 0 ; | |
959 | Float_t meanMassError = 0 ; | |
960 | Float_t meanWidth = 0 ; | |
961 | Float_t meanWidthError = 0 ; | |
962 | Int_t cnt = 0 ; | |
963 | for (Int_t qw = 0 ; qw < nofMassHistograms ; qw++){ | |
53399039 | 964 | if ( fitResultsMean[qw] > invMassMin && fitResultsMean[qw] < invMassMax && fitResultsWidth[qw] > 0 && fitResultsWidth[qw] < 1){ |
29e08f93 | 965 | meanWidth += fitResultsWidth[qw] ; |
966 | meanWidthError += fitResultsWidthErr[qw]; | |
967 | meanMass += fitResultsMean[qw] ; | |
968 | meanMassError += fitResultsMeanErr[qw]; | |
969 | cnt++; | |
970 | } | |
971 | } | |
972 | if (cnt==0) { | |
973 | cout << "Fitting procedure didn't work" << endl; | |
974 | break ; | |
975 | } | |
976 | ||
977 | cout << " Mass : " << meanMass/cnt << " +/- " << meanMassError/cnt << endl ; | |
978 | cout << " Width : " << meanWidth/cnt << " +/- " << meanWidthError/cnt << endl ; | |
979 | cout << " Mean Chi2 of the fits : " << chi2sum/nofMassHistograms << endl ; | |
980 | cout << "********************************************" << endl; | |
981 | cout << " " << endl; | |
982 | ||
983 | ||
984 | cout << " " << endl; | |
985 | cout << "***********************" << endl; | |
53399039 | 986 | cout << " Integrated efficiency (+/- "<< 3*masssigma << " GeV around Resonance mass cut)" << endl ; |
987 | cout << " ResonanceESDAcc/ResonanceMCAcc = " << hptyResonanceESDAcc->GetEntries() << "/" << hptyResonanceMCAcc->GetEntries() << " = " << hptyResonanceESDAcc->GetEntries()/hptyResonanceMCAcc->GetEntries() <<endl; | |
29e08f93 | 988 | cout << "***********************" << endl; |
989 | cout << " " << endl; | |
990 | ||
991 | ||
992 | cout << " " << endl; | |
993 | cout << "***********************" << endl; | |
994 | cout << " Trigger efficiency" << endl ; | |
53399039 | 995 | cout << " Two muons matching = " << hInvMassResonanceTriggerTwoMatch->GetEntries() << "/" << hInvMassResonanceTrigger->GetEntries() << " = " << hInvMassResonanceTriggerTwoMatch->GetEntries()/hInvMassResonanceTrigger->GetEntries() << endl; |
996 | cout << " Single muon matching = " << hInvMassResonanceTriggerOneMatch->GetEntries() << "/" << hInvMassResonanceTrigger->GetEntries() << " = " << hInvMassResonanceTriggerOneMatch->GetEntries()/hInvMassResonanceTrigger->GetEntries() << endl; | |
997 | cout << " No matching = " << hInvMassResonanceTriggerNoMatch->GetEntries() << "/" << hInvMassResonanceTrigger->GetEntries() << " = " << hInvMassResonanceTriggerNoMatch->GetEntries()/hInvMassResonanceTrigger->GetEntries() <<endl; | |
29e08f93 | 998 | cout << "***********************" << endl; |
999 | cout << " " << endl; | |
1000 | ||
1001 | ||
1002 | ||
1003 | ||
1004 | TH2F *hFrame = new TH2F("hFrame", "A hFrame",ptbins,ptmin,ptmax,100,fitResultsWidth[1]-fitResultsWidth[1]*2,fitResultsWidth[1]+fitResultsWidth[1]*2); | |
53399039 | 1005 | hptyResonanceMCAcc->SetLineColor(1); |
1006 | hptyResonanceMCAcc->SetLineStyle(1); | |
1007 | hptyResonanceMCAcc->SetLineWidth(2); | |
29e08f93 | 1008 | |
1009 | ||
53399039 | 1010 | TCanvas *cA = new TCanvas("cA", "Resonance : Width from fit",30,330,700,500); |
29e08f93 | 1011 | cA->Range(-1.69394,-0.648855,15.3928,2.77143); |
1012 | cA->SetBorderSize(2); | |
1013 | cA->SetRightMargin(0.0229885); | |
1014 | cA->SetTopMargin(0.0275424); | |
1015 | cA->SetFrameFillColor(0); | |
1016 | cA->cd(); | |
1017 | ||
1018 | //TH2F *hFrameA = new TH2F("hFrameA", "A hFrameA",ptbins,ptmin,ptmax,100,fitResultsWidth[1]-fitResultsWidth[1]*2,fitResultsWidth[1]+fitResultsWidth[1]*2); | |
1019 | TH2F *hFrameA = new TH2F("hFrameA", "A hFrameA",ptbins,ptmin,ptmax,100,0.,2*fitResultsWidth[1]); | |
1020 | hFrameA->GetYaxis()->SetTitle("Peak Width [GeV/c^{2}]"); | |
1021 | hFrameA->GetYaxis()->SetTitleOffset(1.2); | |
1022 | hFrameA->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
1023 | hFrameA->SetStats(0); | |
1024 | hFrameA->Draw(); | |
1025 | ||
1026 | ||
1027 | TGraphErrors *grWidth = new TGraphErrors(nofMassHistograms,ptbinscenters,fitResultsWidth,errx,fitResultsWidthErr); | |
53399039 | 1028 | grWidth->SetTitle("Resonance Mass Width"); |
29e08f93 | 1029 | grWidth->SetMarkerColor(4); |
1030 | grWidth->SetMarkerStyle(21); | |
1031 | grWidth->Draw("P"); | |
1032 | ||
1033 | ||
1034 | ||
53399039 | 1035 | TCanvas *cB = new TCanvas("cB", "Resonance : Mean from fit",50,350,700,500); |
29e08f93 | 1036 | cB->Range(-1.69394,-0.648855,15.3928,2.77143); |
1037 | cB->SetBorderSize(2); | |
1038 | cB->SetRightMargin(0.0229885); | |
1039 | cB->SetTopMargin(0.0275424); | |
1040 | cB->SetFrameFillColor(0); | |
1041 | cB->cd(); | |
1042 | ||
1043 | //TH2F *hFrameB = new TH2F("hFrameB", "A hFrameB",ptbins,ptmin,ptmax,100,fitResultsMean[1]-fitResultsMean[1]*2,fitResultsMean[1]+fitResultsMean[1]*2); | |
53399039 | 1044 | TH2F *hFrameB = new TH2F("hFrameB", "A hFrameB",ptbins,ptmin,ptmax,100,invMassMin,invMassMax); |
29e08f93 | 1045 | hFrameB->GetYaxis()->SetTitle("Peak Position [GeV/c^{2}]"); |
1046 | hFrameB->GetXaxis()->SetTitle("P_{#perp} [GeV/c]"); | |
1047 | hFrameB->SetStats(0); | |
1048 | hFrameB->Draw(); | |
1049 | ||
1050 | ||
1051 | ||
1052 | ||
1053 | // TGaxis *Mean = new TGaxis(16,fitResultsWidth[1]-fitResultsWidth[1]*2,16,fitResultsWidth[1]+fitResultsWidth[1]*2,7,12,510,"+"); | |
1054 | // Mean->SetLabelOffset(0.050); | |
1055 | // Mean->SetLabelSize(0.04); | |
1056 | // Mean->SetTickSize(0.03); | |
1057 | // Mean->SetGridLength(0); | |
1058 | // Mean->SetTitleOffset(0.1); | |
1059 | // Mean->SetTitleOffset(-1.1); | |
1060 | // Mean->SetTitleSize(0.05); | |
1061 | // Mean->SetTitle("Mean [Gev/c^2]"); | |
1062 | // Mean->Draw(); | |
1063 | ||
1064 | ||
1065 | TGraphErrors *grMean = new TGraphErrors(nofMassHistograms,ptbinscenters,fitResultsMean,errx,fitResultsMeanErr); | |
53399039 | 1066 | grMean->SetTitle("Resonance Mass Mean"); |
29e08f93 | 1067 | grMean->SetMarkerColor(4); |
1068 | grMean->SetMarkerStyle(21); | |
1069 | grMean->Draw("P"); | |
1070 | ||
1071 | cA->Update(); | |
1072 | cB->Update(); | |
1073 | if (SAVE){ | |
53399039 | 1074 | cB->SaveAs("ResonanceMeanVsPt.gif"); |
1075 | cB->SaveAs("ResonanceMeanVsPt.eps"); | |
1076 | cA->SaveAs("ResonanceWidthVsPt.gif"); | |
1077 | cA->SaveAs("ResonanceWidthVsPt.eps"); | |
29e08f93 | 1078 | } |
1079 | ||
1080 | if (WRITE){ | |
1081 | TFile *myFile = new TFile("MUONplotefficiency.root", "RECREATE"); | |
53399039 | 1082 | hptyResonanceMCAcc->Write(); |
1083 | hptyResonanceMC->Write(); | |
1084 | hptResonanceMCAcc ->Write(); | |
1085 | hyResonanceMCAcc->Write(); | |
29e08f93 | 1086 | |
53399039 | 1087 | hptyResonanceESDAcc->Write(); |
1088 | hptyResonanceESD->Write(); | |
1089 | hptResonanceESDAcc ->Write(); | |
1090 | hyResonanceESDAcc->Write(); | |
29e08f93 | 1091 | |
53399039 | 1092 | hEffptyResonance->Write(); |
1093 | hEffyResonance->Write(); | |
1094 | hEffptResonance->Write(); | |
29e08f93 | 1095 | |
1096 | hInvMass->Write(); | |
53399039 | 1097 | hInvMassResonanceTrigger->Write(); |
1098 | hResonanceTriggerNo->Write(); | |
1099 | hInvMassResonanceTriggerOneMatch->Write(); | |
1100 | hInvMassResonanceTriggerTwoMatch->Write(); | |
1101 | hInvMassResonanceTriggerNoMatch->Write(); | |
29e08f93 | 1102 | |
1103 | for (Int_t qw = 0 ; qw < nofMassHistograms ; qw++) | |
1104 | hInvMassInPtBins[qw]->Write(); | |
1105 | ||
1106 | hInvMassAll->Write(); | |
1107 | ||
29e08f93 | 1108 | myFile->Close(); |
1109 | } | |
1110 | ||
1111 | } // macro ends |