]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/Correlations/macros/minijet/analyse_pA.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGCF / Correlations / macros / minijet / analyse_pA.C
CommitLineData
c683985a 1/* $Id: $ */
2//--------------------------------------------------
3//
4// macro to do the final analysis step
5// uses input of analysis class AliAnalysisTaskPhiCorrelation
6//
7// Author : Emilia Leogrande (University of Utrecht)
8//
9//-------------------------------------------------
10
11#include <TChain.h>
12#include <TList.h>
13#include <TTree.h>
14#include <TH1F.h>
15#include <TH2F.h>
16#include <TH3F.h>
17#include <THnSparse.h>
18#include <TProfile.h>
19#include <TCanvas.h>
20#include "TRandom.h"
21#include "TGraphErrors.h"
22#include "TFile.h"
23#include "TF1.h"
24#include "TMath.h"
25#include "TDirectory.h"
26#include "TStyle.h"
27#include "TROOT.h"
28#include "TColor.h"
29
30#include <iostream>
31using namespace std;
32
33void analyseEmy2(Bool_t zyam = kTRUE); // if zyam = kFALSE, fit is used
34Double_t fitFunction(Double_t *x ,Double_t *par); // fit function using constant + 3 gaussians
35Double_t fitFunction2Gaus(Double_t *x ,Double_t *par); // fit function using constant + 2 gaussians
36
37//input file and mixed event removed file
38TFile *fileData=0x0;
39TFile *fileDataEMremoved = 0x0;
40
41const int multclass = 20;
42
43TH1D *fDeltaPhiNch[multclass];
44TH1D *fDeltaEtaNch[multclass];
45TH1D *fSignalDPhi[multclass];
46TH1D *fSignalNSDPhi[multclass];
47TH1D *fSignalASDPhi[multclass];
48TH1D *fRidge1DPhi[multclass];
49TH1D *fRidge2DPhi[multclass];
50TH1D *fRidgeDPhi[multclass];
51TH1D *fSymmRidgeNotScaled[multclass];
52TH1D *fSymmRidge[multclass];
53TH1D *fFinal1DPhi[multclass];
54TH1D *fFinalDPhi[multclass];
55
56TString flag = "R";
57TF1 *fTotal2Gaus[multclass]; // fit with 2 gaussians + const
58TF1 *fTotal[multclass]; // fit with 3 gaussians + const
59
60//properties of histogram
61const int bins = 72; //
62Double_t binWidth=2*TMath::Pi()/bins;
63
64const int binsDeta = 48;
65
66
67Double_t max_bin_for_etagap = 1.2;
68Double_t min_bin_for_etagap = -1.2;
69Double_t max_eta = 1.8;
70Double_t min_eta = -1.8;
71
72//________________________________________________________________________________________________________________
73//
74Double_t fitFunction(Double_t *x ,Double_t *par)
75{
76 // fit function for 3 gaus + constant
77
78 // parameters for Gaussian
79 Double_t A1 = par[0];
80 Double_t sigma1 = par[1];
81 Double_t A2 = par[2];
82 Double_t sigma2 = par[3];
83 Double_t A3 = par[4];
84 Double_t sigma3 = par[5];
85 Double_t integral = par[6];
86
87 Double_t constante = (integral-
88 TMath::Sqrt(TMath::Pi()*2)/ binWidth*
89 (A1 * sigma1 + A2 * sigma2 + A3*sigma3))/bins;
90 Double_t q = x[0];
91
92 //fit value
93 Double_t fitval = constante +
94 (q>-0.5*TMath::Pi()&&q<0.5*TMath::Pi())*(
95 A1 * exp(- q * q / (2 * sigma1 *sigma1)) +
96 A1 * exp(-((q - TMath::TwoPi())) * ((q - TMath::TwoPi())) / ( 2 * sigma1 * sigma1))
97 )
98 +
99 (q>-0.2*TMath::Pi()&&q<0.2*TMath::Pi())*(
100 A2 * exp(- q * q / (2 * sigma2 *sigma2)) +
101 A2 * exp(-((q - TMath::TwoPi())) * ((q - TMath::TwoPi())) / ( 2 * sigma2 * sigma2))
102 )
103 +
104 (q>0.5*TMath::Pi()&&q<1.5*TMath::Pi())*(
105 A3 * exp(-((q - TMath::Pi())) * ((q - TMath::Pi())) / ( 2 * sigma3 * sigma3)) +
106 A3 * exp(-((q + TMath::Pi())) * ((q + TMath::Pi())) / (2 * sigma3 * sigma3))
107 );
108 return fitval;
109}
110
111//________________________________________________________________________________________________________________
112//
113Double_t fitFunction2Gaus(Double_t *x ,Double_t *par)
114{
115 // fit function for 2 gaus + constant
116
117 // parameters for Gaussian
118 Double_t A1 = par[0];
119 Double_t sigma1 = par[1];
120 Double_t A3 = par[2];
121 Double_t sigma3 = par[3];
122 Double_t integral = par[4];
123
124 Double_t constante = (integral -
125 TMath::Sqrt(TMath::Pi()*2)/ binWidth*
126 (A1 * sigma1 + A3*sigma3))/bins;
127 Double_t q = x[0];
128
129 //fit value
130 Double_t fitval = constante +
131 (q>-0.5*TMath::Pi()&&q<0.5*TMath::Pi())*(
132 A1 * exp(- q * q / (2 * sigma1 *sigma1)) +
133 A1 * exp(-((q - TMath::TwoPi())) * ((q - TMath::TwoPi())) / ( 2 * sigma1 * sigma1))
134 )
135 +
136 (q>0.5*TMath::Pi()&&q<1.5*TMath::Pi())*(
137 A3 * exp(-((q - TMath::Pi())) * ((q - TMath::Pi())) / ( 2 * sigma3 * sigma3)) +
138 A3 * exp(-((q + TMath::Pi())) * ((q + TMath::Pi())) / (2 * sigma3 * sigma3))
139 );
140 return fitval;
141}
142
143//_______________________________________________________________________________________________________________
144//
145Double_t fline(Double_t *x, Double_t *par){
146
147 if(x[0]>-1.8 && x[0]<=0){
148 return par[0]+par[1]*x[0];
149 }
150 else if(x[0]>0 && x[0]<1.8){
151 return par[2]+par[3]*x[0];
152 }
153 else
154 return 0;
155}
156
157
158//________________________________________________________________________________________________________________
159//
160void analyseEmy2(Bool_t zyam){
161
162
163 // plot style
164 gStyle->SetOptStat(0);
165 const Int_t NRGBs = 5;
166 const Int_t NCont = 500;
167 Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
168 Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
169 Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
170 Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
171 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
172 gStyle->SetNumberContours(NCont);
173
174 //style
175 gROOT->SetStyle("Plain");
176 gStyle->SetOptStat(0);
177 gStyle->SetPalette(1);
178
179 //-------------- TRIGGERS AND EVENTS
180
181 TH2D *dphideta[multclass];
182 TH1D * trigger = 0x0;
183 TH1D * event = 0x0;
184
185 fileData = TFile::Open("dphi_corr.root");
186 trigger = (TH1D*)fileData->Get("triggers_0");
187 event = (TH1D*)fileData->Get("events");
188
189 // get average trigger particles per event
190 TProfile *p0 = (TProfile*)trigger->Clone();
191 TProfile *p1 = (TProfile*)event->Clone();
192 p0->Sumw2();
193 p1->Sumw2();
194 p0->Divide(p0,p1,1,1,"B");
195
196 // copy triggers and events in the new dphi_corr with the Mixed Event removed
197 TH1D *triggerCopy = 0x0;
198 TH1D *eventCopy = 0x0;
199
200 triggerCopy = (TH1D*)trigger->Clone();
201 eventCopy = (TH1D*)event->Clone();
202
203 fileDataEMremoved = TFile::Open("dphi_corr_MEremoved.root","RECREATE");
204 triggerCopy->SetName("triggers_0");
205 triggerCopy->Write();
206 eventCopy->SetName("events");
207 eventCopy->Write();
208 fileDataEMremoved->Close();
209
210
211 //-------------- MIXED EVENT REMOVAL: restores the right number of particles in the detector acceptance but keeps the detector azimuthal unefficiencies corrections and cures the dip in (0,0) from two-trak cuts
212 // Removing the event mixing: S/M (from dphi_corr) * M (from the triangle)
213
214 Double_t triangle_factor[binsDeta]={0};
215
216 TH2D *s_over_m[multclass];
217 TH1D *s_m_deta[multclass];
218 TH2D *s_over_m_x_m[multclass];
219
220 for(Int_t i=0;i<multclass;i++){
221 s_over_m[i] = (TH2D*)fileData->Get(Form("dphi_0_0_%d",i));
222 s_m_deta[i] = (TH1D*)s_over_m[i]->ProjectionY()->Clone();
223 s_over_m_x_m[i] = (TH2D*)s_over_m[i]->Clone();
224 s_over_m_x_m[i]->Reset();
225 }
226
227
228 TF1 *f2 = new TF1("f2",fline,min_eta,max_eta,4);
229
230 f2->FixParameter(0,1);
231 f2->FixParameter(1,1/max_eta);
232 f2->FixParameter(2,1);
233 f2->FixParameter(3,-1/max_eta);
234
235 for(Int_t i=0;i<binsDeta;i++){
236
237 triangle_factor[i] = f2->Eval(s_m_deta[0]->GetBinCenter(i+1));
238
239 }
240
241
242
243 //--scale each deta bin of the old TH2 with the triangle_factor[deta]
244
245 for(Int_t i=0;i<multclass;i++){
246 for(Int_t j=0;j<binsDeta;j++){
247 for(Int_t k=0;k<bins;k++){
248 s_over_m_x_m[i] -> SetBinContent(k+1,j+1,(s_over_m[i]->GetBinContent(k+1,j+1))*triangle_factor[j]);
249 s_over_m_x_m[i]->SetBinError(k+1,j+1,(s_over_m[i]->GetBinError(k+1,j+1))*triangle_factor[j]);
250 }
251 }
252 }
253
254 fileDataEMremoved = TFile::Open("dphi_corr_MEremoved.root","UPDATE");
255
256 for(Int_t i=0;i<multclass;i++){
257
258 s_over_m_x_m[i]->SetName(Form("dphiNoMixed_%d",i));
259 s_over_m_x_m[i]->Write();
260
261 }
262
263
264
265 //-------------- DOUBLE RIDGE SUBTRACTION: gets rid of no-jet related components (v3 is still kept => effect added to the systematics)
266
267 // the ridge, estimated via an etagap, has to be scaled since it sits on the triangle
268 Double_t scale_for_ridge_NS = 0, scale_for_ridge_AS = 0;
269
270
271 scale_for_ridge_NS = f2->Integral(min_bin_for_etagap,max_bin_for_etagap)/(f2->Integral(min_eta,min_bin_for_etagap)+f2->Integral(max_bin_for_etagap,max_eta)); //there is etagap in the NS
272 cout<<"scaling NS:"<<scale_for_ridge_NS<<endl;
273
274 scale_for_ridge_AS = f2->Integral(min_eta,max_eta)/(f2->Integral(min_eta,min_bin_for_etagap)+f2->Integral(max_bin_for_etagap,max_eta)); // there is no etagap in the AS
275 cout<<"scaling AS:"<<scale_for_ridge_AS<<endl;
276
277 // Double ridge subtraction
278
279 TCanvas *c = new TCanvas();
280 c->Divide(5,4);
281
282 for(Int_t i=0;i<multclass;i++){
283 c->cd(i+1);
284
285
286 dphideta[i] = (TH2D*)fileDataEMremoved->Get(Form("dphiNoMixed_%d",i));
287
288
289 // phi and eta projections
290 fDeltaPhiNch[i] = (TH1D*)dphideta[i]->ProjectionX()->Clone();
291 if(!zyam)
292 fDeltaPhiNch[i]->Scale(binWidth); //gaussians include the binwidth, so when using the fit, the histograms must be scaled first
293 fDeltaPhiNch[i]->Draw();
294
295 fDeltaEtaNch[i] = (TH1D*)dphideta[i]->ProjectionY()->Clone();
296
297 // signal NS: |DEta|<max_bin_for_etagap; signal AS: |DEta|<max_eta
298 fSignalNSDPhi[i] = (TH1D*)dphideta[i]->ProjectionX(Form("|DEta|<%f",max_bin_for_etagap),fDeltaEtaNch[i]->FindBin(min_bin_for_etagap+0.0001),fDeltaEtaNch[i]->FindBin(max_bin_for_etagap-0.0001))->Clone();
299 fSignalASDPhi[i] = (TH1D*)dphideta[i]->ProjectionX(Form("|DEta|<%f",max_eta))->Clone();
300
301 fSignalDPhi[i] = (TH1D*)fSignalASDPhi[i]->Clone();
302 fSignalDPhi[i]->Reset();
303 fSignalDPhi[i]->Sumw2();
304
305 for(Int_t k=0;k<bins/2;k++){
306 fSignalDPhi[i]->SetBinContent(k+1,fSignalNSDPhi[i]->GetBinContent(k+1));
307 fSignalDPhi[i]->SetBinError(k+1, fSignalNSDPhi[i]->GetBinError(k+1));
308 }
309 for(Int_t k=bins/2;k<bins;k++){
310 fSignalDPhi[i]->SetBinContent(k+1,fSignalASDPhi[i]->GetBinContent(k+1));
311 fSignalDPhi[i]->SetBinError(k+1, fSignalASDPhi[i]->GetBinError(k+1));
312 }
313 if(!zyam)
314 fSignalDPhi[i]->Scale(binWidth);
315
316 // ridge1 DEta<min_bin_for_etagap
317 fRidge1DPhi[i] = (TH1D*)dphideta[i]->ProjectionX(Form("DEta<%f",min_bin_for_etagap),1,fDeltaEtaNch[i]->FindBin(min_bin_for_etagap-0.0001))->Clone();
318 if(!zyam)
319 fRidge1DPhi[i]->Scale(binWidth);
320 fRidge1DPhi[i]->SetMarkerColor(kRed);
321
322 // ridge2 DEta>max_bin_for_etagap
323 fRidge2DPhi[i] = (TH1D*)dphideta[i]->ProjectionX(Form("DEta>%f",max_bin_for_etagap),fDeltaEtaNch[i]->FindBin(max_bin_for_etagap+0.0001),fDeltaEtaNch[i]->GetNbinsX())->Clone();
324 if(!zyam)
325 fRidge2DPhi[i]->Scale(binWidth);
326 fRidge2DPhi[i]->SetMarkerColor(kBlue);
327
328 // ridge = ridge1 + ridge2
329 fRidgeDPhi[i] = (TH1D*)fRidge1DPhi[i]->Clone("fRidge");
330 fRidgeDPhi[i]->Reset();
331 fRidgeDPhi[i]->Sumw2();
332 fRidgeDPhi[i]->Add(fRidge1DPhi[i],fRidge2DPhi[i],1,1);
333 //fRidgeDPhi[i]->Scale(scale_for_ridge);
334
335 // symmetrize NS ridge in the AS
336 fSymmRidgeNotScaled[i] = (TH1D*)fRidgeDPhi[i]->Clone("fSymmRidgeNotScaled");
337
338 for(Int_t k=fSymmRidgeNotScaled[i]->GetNbinsX()/2+1;k<=fSymmRidgeNotScaled[i]->GetNbinsX();k++){
339
340 fSymmRidgeNotScaled[i]->SetBinContent(k,fSymmRidgeNotScaled[i]->GetBinContent(fSymmRidgeNotScaled[i]->GetNbinsX()+1-k));
341
342 }
343
344 // scale the symmetrized ridge according to NS or AS
345 fSymmRidge[i] = (TH1D*)fSymmRidgeNotScaled[i]->Clone("fSymmRidge");
346
347 for(Int_t k=0;k<bins/2;k++){
348 fSymmRidge[i]->SetBinContent(k+1,(fSymmRidgeNotScaled[i]->GetBinContent(k+1))*scale_for_ridge_NS);
349 }
350 for(Int_t k=bins/2;k<bins;k++){
351 fSymmRidge[i]->SetBinContent(k+1,(fSymmRidgeNotScaled[i]->GetBinContent(k+1))*scale_for_ridge_AS);
352 }
353
354
355 // signal - symmetric ridge
356
357 if(zyam){
358 fFinal1DPhi[i] = new TH1D(Form("fFinal1DPhi[%d]",i),Form("fFinal1DPhi[%d]",i),bins,-0.5*TMath::Pi(),1.5*TMath::Pi());
359 fFinal1DPhi[i]->Add(fSignalDPhi[i],fSymmRidge[i],1,-1);
360 fFinal1DPhi[i]->Sumw2();
361 fFinalDPhi[i] = (TH1D*)fFinal1DPhi[i]->Clone("fFinal"); // zyam: average between the two min values => sum first half of NS in the second half and second half of AS in the first half, so zyam = min/2
362 fFinalDPhi[i]->Reset();
363 fFinalDPhi[i]->Sumw2();
364
365 for(Int_t k=1;k<=bins/4;k++){
366 fFinalDPhi[i]->SetBinContent(k,0.);
367 fFinalDPhi[i]->SetBinContent(k+bins/4,fFinal1DPhi[i]->GetBinContent(k+bins/4)+fFinal1DPhi[i]->GetBinContent(bins/4+1-k));
368 fFinalDPhi[i]->SetBinError(k+bins/4,TMath::Sqrt(pow(fFinal1DPhi[i]->GetBinError(k+bins/4),2)+pow(fFinal1DPhi[i]->GetBinError(bins/4+1-k),2)));
369 fFinalDPhi[i]->SetBinContent(k+bins/2,fFinal1DPhi[i]->GetBinContent(k+bins/2)+fFinal1DPhi[i]->GetBinContent(bins+1-k));
370 fFinalDPhi[i]->SetBinError(k+bins/2,TMath::Sqrt(pow(fFinal1DPhi[i]->GetBinError(k+bins/2),2)+pow(fFinal1DPhi[i]->GetBinError(bins+1-k),2)));
371 fFinalDPhi[i]->SetBinContent(k+bins/4*3,0.);
372
373 }
374 }
375
376 else{
377
378 fFinalDPhi[i] = (TH1D*)fSignalDPhi[i]->Clone();
379 fFinalDPhi[i]->Reset();
380 fFinalDPhi[i]->Sumw2();
381 fFinalDPhi[i]->Add(fSignalDPhi[i],fSymmRidge[i],1,-1);
382 }
383
384 }
385
386 // store the pair yields in a file (the yields are *not* normalized to the Ntriggers)
387
388 TFile* file_yields = 0x0;
389 if(zyam)
390 file_yields = TFile::Open("PairYields_zyam.root","RECREATE");
391 else
392 file_yields = TFile::Open("PairYields_fit.root","RECREATE");
393
394
395 for(Int_t i=0;i<multclass;i++){
396 fDeltaEtaNch[i]->SetName(Form("DeltaEta_0_0_%d",i));
397 fDeltaEtaNch[i]->Write();
398 fDeltaPhiNch[i]->SetName(Form("Correlation bin %d in dphi",i));
399 fDeltaPhiNch[i]->Write();
400 fSignalDPhi[i]->SetName(Form("Signal_0_0_%d",i));
401 fSignalDPhi[i]->Write();
402 fRidgeDPhi[i]->SetName(Form("Ridge_0_0_%d",i));
403 fRidgeDPhi[i]->Write();
404 fSymmRidgeNotScaled[i]->SetName(Form("Symmetric_Ridge_NotScaled_0_0_%d",i));
405 fSymmRidgeNotScaled[i]->Write();
406 fSymmRidge[i]->SetName(Form("Symmetric_Ridge_0_0_%d",i));
407 fSymmRidge[i]->Write();
408 fFinalDPhi[i]->SetName(Form("Pure_Signal_0_0_%d",i));
409 fFinalDPhi[i]->Write();
410 }
411 file_yields->Close();
412
413 //-------------- CORRELATION OBSERVABLES: per-trigger yields, triggers and uncorrelated seeds
414
415 Float_t baseline[multclass]={0};
416
417 TGraphErrors *fNearSideIntegral = new TGraphErrors();
418 fNearSideIntegral->SetName("fNearSideIntegral");
419 fNearSideIntegral->SetMarkerColor(kGreen+2);
420 fNearSideIntegral->SetLineColor(kGreen+2);
421 fNearSideIntegral->SetLineWidth(1);
422 fNearSideIntegral->SetMarkerStyle(4);
423
424 TGraphErrors *fAwaySideIntegral = new TGraphErrors();
425 fAwaySideIntegral->SetName("fAwaySideIntegral");
426 fAwaySideIntegral->SetMarkerColor(kBlue);
427 fAwaySideIntegral->SetLineColor(kBlue);
428 fAwaySideIntegral->SetLineWidth(1);
429 fAwaySideIntegral->SetMarkerStyle(4);
430
431 TGraphErrors *fBothSideIntegral = new TGraphErrors();
432 fBothSideIntegral->SetName("fBothSideIntegral");
433 fBothSideIntegral->SetMarkerColor(kMagenta);
434 fBothSideIntegral->SetLineColor(kMagenta);
435 fBothSideIntegral->SetLineWidth(1);
436 fBothSideIntegral->SetMarkerStyle(4);
437
438
439 TGraphErrors *fNjets = new TGraphErrors();
440 fNjets->SetName("fNjets");
441 fNjets->SetMarkerColor(kCyan+2);
442 fNjets->SetLineColor(kCyan+2);
443 fNjets->SetLineWidth(1);
444 fNjets->SetMarkerStyle(4);
445
446 TGraphErrors *fTriggerAverage = new TGraphErrors();
447 fTriggerAverage->SetName("fTriggerAverage");
448 fTriggerAverage->SetMarkerColor(kBlack);
449 fTriggerAverage->SetLineColor(kBlack);
450 fTriggerAverage->SetLineWidth(1);
451 fTriggerAverage->SetMarkerStyle(4);
452
453 Int_t points=0;
454 Double_t minbin[multclass] = {0};
455
456 // extract information out of dphi histograms
457 TCanvas * cYields= new TCanvas("cYields", "cYields", 150, 150, 820, 620);
458 cYields->Divide(5,4);
459
460 for(Int_t i=0;i<multclass;i++){
461 cYields->cd(i+1);
462
463
464 if(zyam) {
465
466 if(fFinalDPhi[i]->Integral()>0){
467 fFinalDPhi[i]->GetXaxis()->SetRange(bins/4+1,bins/4*3);
468 baseline[i]=fFinalDPhi[i]->GetMinimum()/2;
469 minbin[i] = fFinalDPhi[i]->GetMinimumBin();
470 fFinalDPhi[i]->GetXaxis()->UnZoom();
471
472 for(Int_t k=0;k<bins;k++){
473 if(fFinalDPhi[i]->GetBinContent(k+1)!=0)
474 fFinalDPhi[i]->SetBinContent(k+1,fFinalDPhi[i]->GetBinContent(k+1)-baseline[i]);
475 else
476 fFinalDPhi[i]->SetBinContent(k+1,0.);
477 }
478
479 fFinalDPhi[i]->DrawClone("");
480
481 fFinalDPhi[i]->SetTitle(Form("0.7<p_{T,trig}<5.0 - 0.7<p_{T,assoc}<5.0 - %d-%d %",i*5,(i+1)*5));
482 fFinalDPhi[i]->SetTitle("1/N_{trig} dN_{assoc}/d#Delta#varphi (rad^{-1})");
483 //-
484 Double_t errorNS = 0;
485 Double_t nearSideResult = (fFinalDPhi[i]->IntegralAndError(0,minbin[i],errorNS,"width"))/trigger->GetBinContent(i+1);
486 Double_t nearSideError = errorNS/trigger->GetBinContent(i+1);
487 fNearSideIntegral->SetPoint(points,i, nearSideResult);
488 fNearSideIntegral->SetPointError(points,0.5,errorNS/trigger->GetBinContent(i+1));
489 //-
490
491 //--
492 Double_t errorAS = 0;
493 Double_t awaySideResult = (fFinalDPhi[i]->IntegralAndError(minbin[i],bins,errorAS,"width"))/trigger->GetBinContent(i+1);
494 Double_t awaySideError = errorAS/trigger->GetBinContent(i+1);
495 fAwaySideIntegral->SetPoint(points,i, awaySideResult );
496 fAwaySideIntegral->SetPointError(points,0.5, errorAS/trigger->GetBinContent(i+1));
497 //--
498
499 //---
500 Double_t bothSideResult = nearSideResult + awaySideResult;
501 Double_t bothSideError = bothSideResult * TMath::Sqrt(pow(errorNS,2)+pow(errorAS,2))/trigger->GetBinContent(i+1);
502 fBothSideIntegral->SetPoint(points,i, bothSideResult );
503 fBothSideIntegral->SetPointError(points,0.5, bothSideError );
504 //---
505
506
507
508 }
509 else{
510 fNearSideIntegral->SetPoint(points,i, 0);
511 fAwaySideIntegral->SetPoint(points,i, 0);
512 fBothSideIntegral->SetPoint(points,i,0);
513 }
514 Double_t p0BinContent=p0->GetBinContent(i+1);
515 Double_t p0BinError=p0->GetBinError(i+1);
516
517 //--------
518 Double_t njets = p0BinContent/(1+bothSideResult);
519 Double_t njetsError = njets*TMath::Sqrt(bothSideError*bothSideError/(1+bothSideResult)/(1+bothSideResult)+p0BinError*p0BinError/p0BinContent/p0BinContent);
520 fNjets->SetPoint(points,i, njets );
521 fNjets->SetPointError(points,0.5,njetsError );
522
523 //-------
524
525 fTriggerAverage->SetPoint(points,i, p0BinContent);
526 fTriggerAverage->SetPointError(points,0.5, p0BinError);
527
528 }
529
530 else if (!zyam){
531
532 if(fFinalDPhi[i]->Integral()>0){
533
534 //first fit function: 2 gauss + const
535 fTotal2Gaus[i] = new TF1(Form("gaus3and2_%d",i), fitFunction2Gaus , -0.5*TMath::Pi(), 1.5*TMath::Pi(), 5);
536 fTotal2Gaus[i]->SetName(Form("gaus3_%d",i));
537 fTotal2Gaus[i]->SetParNames ("A1","sigma1","A3", "sigma3");
538 fTotal2Gaus[i]->SetLineColor(kRed);
539 fTotal2Gaus[i]->SetLineWidth(2);
540
541 baseline[i]=fFinalDPhi[i]->GetMinimum();
542 Double_t integr_for_const_2 = fFinalDPhi[i]->Integral();
543
544 fTotal2Gaus[i]->FixParameter(4,integr_for_const_2);
545 fTotal2Gaus[i]->SetParameters( fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(0)) - baseline[i] , 0.6 , fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(TMath::Pi()))-baseline[i] , 0.6);
546
547 fTotal2Gaus[i]->SetParLimits(0, 0, (fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(0))-baseline[i])*2);
548 fTotal2Gaus[i]->SetParLimits(1, 0.01, 10);
549 fTotal2Gaus[i]->SetParLimits(2, 0, (fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(TMath::Pi()))-baseline[i])*2);
550 fTotal2Gaus[i]->SetParLimits(3, 0.01, 10);
551
552 fTotal2Gaus[i]->SetLineColor(kRed);
553 fTotal2Gaus[i]->SetLineWidth(2);
554
555 fFinalDPhi[i]->Fit(fTotal2Gaus[i],flag);
556 fFinalDPhi[i]->SetMinimum(0);
557 fFinalDPhi[i]->DrawClone("");
558 fTotal2Gaus[i] ->DrawClone("same");
559
560 Double_t A11 = fTotal2Gaus[i]->GetParameter(0);
561 Double_t sigma11 = fTotal2Gaus[i]->GetParameter(1);
562 Double_t A31 = fTotal2Gaus[i]->GetParameter(2);
563 Double_t sigma31 = fTotal2Gaus[i]->GetParameter(3);
564
565 Double_t a1e1 = fTotal2Gaus[i]->GetParError(0);
566 Double_t s1e1 = fTotal2Gaus[i]->GetParError(1);
567 Double_t a3e1 = fTotal2Gaus[i]->GetParError(2);
568 Double_t s3e1 = fTotal2Gaus[i]->GetParError(3);
569
570
571 Double_t T11 = A11*sigma11;
572 Double_t T31 = A31*sigma31;
573 Double_t t11 = T11*TMath::Sqrt(a1e1*a1e1/A11/A11 + s1e1*s1e1/sigma11/sigma11);
574 Double_t t31 = T31*TMath::Sqrt(a3e1*a3e1/A31/A31 + s3e1*s3e1/sigma31/sigma31);
575
576
577 //second fit: 3 gauss + const
578 fTotal[i] = new TF1(Form("gaus3_%d",i), fitFunction , -0.5*TMath::Pi(), 1.5*TMath::Pi(), 7);
579 fTotal[i]->SetName(Form("gaus3_%d",i));
580 fTotal[i]->SetParNames ("A1","sigma1","A2","sigma2", "A3", "sigma3","integral");
581 fTotal[i]->SetLineColor(kRed);
582 fTotal[i]->SetLineWidth(2);
583
584 Double_t integr_for_const = fFinalDPhi[i]->Integral();
585
586
587 fTotal[i]->FixParameter(0,A11);
588 fTotal[i]->FixParameter(1,sigma11*1.2);
589 fTotal[i]->FixParameter(2,A11);
590 fTotal[i]->FixParameter(3,sigma11*0.7);
591 fTotal[i]->FixParameter(4,A31);
592 fTotal[i]->FixParameter(5,sigma31);
593 fTotal[i]->FixParameter(6,integr_for_const);
594
595 fTotal[i]->SetParLimits(0, 0, (fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(0))-baseline[i])*2);
596 fTotal[i]->SetParLimits(1, 0.3, 10);
597 fTotal[i]->SetParLimits(2, 0, (fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(0))-baseline[i])*2);
598 fTotal[i]->SetParLimits(3, 0.12, 0.4);
599 fTotal[i]->SetParLimits(4, 0, (fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->
600 GetXaxis()->FindFixBin(TMath::Pi()))-baseline[i])*2);
601 fTotal[i]->SetParLimits(5, 0.01, 10);
602
603 fTotal[i]->SetLineColor(kRed);
604 fTotal[i]->SetLineWidth(2);
605
606
607 fFinalDPhi[i]->Fit(fTotal[i],flag);
608 fFinalDPhi[i]->SetMinimum(0);
609 fFinalDPhi[i]->DrawClone("");
610 fFinalDPhi[i]->SetTitle(Form("0.7<p_{T,trig}<5.0 - 0.7<p_{T,assoc}<5.0 - %d-%d %",i*5,(i+1)*5));
611 fFinalDPhi[i]->SetTitle("1/N_{trig} dN_{assoc}/d#Delta#varphi (rad^{-1})");
612 fTotal[i]->DrawClone("same");
613
614 Double_t A1 = fTotal[i]->GetParameter(0);
615 Double_t sigma1 = fTotal[i]->GetParameter(1);
616 Double_t A2 = fTotal[i]->GetParameter(2);
617 Double_t sigma2 = fTotal[i]->GetParameter(3);
618 Double_t A3 = fTotal[i]->GetParameter(4);
619 Double_t sigma3 = fTotal[i]->GetParameter(5);
620
621
622 //define each gaussian and constant to be drawn with different colors on top of each other
623
624 TF1 * fConstant = new TF1("konst", "pol0(0)",-0.5*TMath::Pi(), 1.5*TMath::Pi());
625 fConstant->SetParameter(0,(integr_for_const - TMath::Sqrt(TMath::Pi()*2)/binWidth*(A1*sigma1+A2*sigma2+A3*sigma3))/bins);
626 fConstant->SetLineColor(kBlue);
627 fConstant->Draw("same");
628
629 //gaus 1 NS
630 TF1 * fGaussian1 = new TF1("fGaussian1", "[0]*exp(-x*x/(2*[1]*[1])) +[0] * exp(-(x-TMath::TwoPi())*(x-TMath::TwoPi())/(2*[1]*[1]))",-0.5*TMath::Pi(), 1.5*TMath::Pi());
631 fGaussian1->SetParameters(fTotal[i]->GetParameter(0),fTotal[i]->GetParameter(1));
632 fGaussian1->SetLineColor(kMagenta);
633 fGaussian1->SetLineStyle(1);
634 fGaussian1->Draw("same");
635
636 //gaus 2 NS
637 TF1 * fGaussian2 = new TF1("fGaussian2", "[0]*exp(-x*x/(2*[1]*[1])) +[0] * exp(-(x-TMath::TwoPi())*(x-TMath::TwoPi())/(2*[1]*[1]))",-0.5*TMath::Pi(), 1.5*TMath::Pi());
638 fGaussian2->SetLineColor(kGreen+2);
639 fGaussian2->SetParameters(fTotal[i]->GetParameter(2),fTotal[i]->GetParameter(3));
640 fGaussian2->Draw("same");
641
642 //gaus 3 AS
643 TF1 * fGaussian3 = new TF1("fGaussian3", "[0] * exp(-((x-TMath::Pi()))*((x-TMath::Pi()))/(2*[1]*[1]))+[0] * exp(-((x+TMath::Pi()))*((x+TMath::Pi()))/(2*[1]*[1]))",-0.5*TMath::Pi(), 1.5*TMath::Pi());
644 fGaussian3->SetLineColor(kCyan);
645 fGaussian3->SetParameters(fTotal[i]->GetParameter(4), fTotal[i]->GetParameter(5));
646 fGaussian3->Draw("same");
647
648
649 Double_t a1e = fTotal[i]->GetParError(0);
650 Double_t s1e = fTotal[i]->GetParError(1);
651 Double_t a2e = fTotal[i]->GetParError(2);
652 Double_t s2e = fTotal[i]->GetParError(3);
653 Double_t a3e = fTotal[i]->GetParError(4);
654 Double_t s3e = fTotal[i]->GetParError(5);
655
656 Double_t T1 = A1*sigma1;
657 Double_t T2 = A2*sigma2;
658 Double_t T3 = A3*sigma3;
659 Double_t t1 = T1*TMath::Sqrt(a1e*a1e/A1/A1 + s1e*s1e/sigma1/sigma1);
660 Double_t t2 = T2*TMath::Sqrt(a2e*a2e/A2/A2 + s2e*s2e/sigma2/sigma2);
661 Double_t t3 = T3*TMath::Sqrt(a3e*a3e/A3/A3 + s3e*s3e/sigma3/sigma3);
662
663 //-
664 Double_t nearSideResult = TMath::Sqrt(TMath::Pi()*2)/ binWidth* (A1 * sigma1 + A2 * sigma2)/trigger->GetBinContent(i+1);
665 Double_t nearSideError = nearSideResult * TMath::Sqrt((t1*t1 + t2*t2)/(T1+T2)/(T1+T2)+ 1./trigger->GetBinContent(i+1));
666 fNearSideIntegral->SetPoint(points,i, nearSideResult);
667 fNearSideIntegral->SetPointError(points,0.5,nearSideError);
668
669 //-
670
671 //--
672 Double_t awaySideResult = TMath::Sqrt(TMath::Pi()*2)/ binWidth*
673 (A3 * sigma3)/trigger->GetBinContent(i+1);
674 Double_t awaySideError = awaySideResult*TMath::Sqrt(a3e*a3e/A3/A3 + s3e*s3e/sigma3/sigma3 + 1/trigger->GetBinContent(i+1));
675 fAwaySideIntegral->SetPoint(points,i, awaySideResult );
676 fAwaySideIntegral->SetPointError(points,0.5, awaySideError );
677 //--
678
679 //---
680 bothSideResult = TMath::Sqrt(TMath::Pi()*2)/ binWidth* (A1 * sigma1 + A2 * sigma2 + A3 * sigma3 )/trigger->GetBinContent(i+1);
681 bothSideError = nearSideResult * TMath::Sqrt((t1*t1 + t2*t2 + t3*t3)/(T1+T2+T3)/(T1+T2+T3)+ 1./trigger->GetBinContent(i+1));
682 fBothSideIntegral->SetPoint(points,i, bothSideResult );
683 fBothSideIntegral->SetPointError(points,0.5, bothSideError );
684 //---
685
686 }
687 else{
688
689 fNearSideIntegral->SetPoint(points,i, 0);
690 fAwaySideIntegral->SetPoint(points,i, 0);
691 fBothSideIntegral->SetPoint(points,i,0);
692
693 }
694 Double_t p0BinContent=p0->GetBinContent(i+1);
695 Double_t p0BinError=p0->GetBinError(i+1);
696
697 //--------
698 Double_t njets = p0BinContent/(1+bothSideResult);
699 Double_t njetsError = njets*TMath::Sqrt(bothSideError*bothSideError/(1+bothSideResult)/(1+bothSideResult) + p0BinError*p0BinError/p0BinContent/p0BinContent);
700 fNjets->SetPoint(points,i, njets );
701 fNjets->SetPointError(points,0.5,njetsError );
702 //-------
703
704 fTriggerAverage->SetPoint(points,i, p0BinContent);
705 fTriggerAverage->SetPointError(points,0.5, p0BinError);
706
707
708 }
709 points++;
710 }
711
712
713 TFile* file = 0x0;
714 if(zyam)
715 file = TFile::Open("njet_zyam.root","RECREATE");
716 else
717 file = TFile::Open("njet_fit.root","RECREATE");
718
719 fNearSideIntegral->Write();
720 fAwaySideIntegral->Write();
721 fBothSideIntegral->Write();
722 fNjets->Write();
723 fTriggerAverage->Write();
724
725 file->Close();
726
727
728
729}
730
731