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