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