In AliMUONTriggerQADataMakerRec:
[u/mrichter/AliRoot.git] / MUON / MUONplotefficiency.C
CommitLineData
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
49Double_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
56Double_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}
62Double_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
84Int_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