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