StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
exampleFits.C
1 // Note: I cut and paste parts of this file into a root session.
2 // I consider this and associated files as templates which can be modified
3 // to suit particular analyses/fitting.
4 
5 // I normally run the root session in a directory which has soft links to the root files.
6 // Output written to a subdirectory fitResults.
7 
8 //Choose dataset to fit.
9 char *dataset = "AuAu200GeV";
10 char *dataset = "AuAu62GeV";
11 char *dataset = "CuCu200GeV";
12 char *dataset = "CuCu62GeV";
13 
14 // Choose cutBin and charge to fit as well as setting number of parameters.
15 int ibin = 0;
16 int icharge = 3;
17 int nParams = 13;
18 
19 // Open appropriate pair of root files.
20 // norm scales from the dN/d\eta observed in the EStruct analysis and the corrected
21 // dN/d\eta from the bulk working group.
22 // fPileup is the fraction of pileup not removed by my pileup cut.
23 // In principle this can be centrality dependent, but I don't know how to measure that.
24 {
25  if (!strcmp(dataset,"AuAu200GeV")) {
26  int nCent = 11;
27  TFile *data1 = new TFile("AuAu200GeV_11c_pair.root");
28  TFile *data2 = new TFile("AuAu200GeV_11c_pair_noPileupCuts.root");
29  const char* centName[] = {"90-100", "80-90", "70-80", "60-70", "50-60", "40-50", "30-40", "20-30", "10-20", "5-10", "0-5"};
30  double norm[] = {1.235, 1.177, 1.157, 1.167, 1.178, 1.193, 1.230, 1.292, 1.366, 1.379, 1.451};
31  // Normalization implied by taking difference between my noPileup analysis and Michael's analysis
32  // and minimizing structure. 10% sort of differences.
33  double normN[] = {1.23, 1.177, 1.11, 1.13, 1.122, 1.142, 1.18, 1.258, 1.364, 1.5, 1.46};
34  double fPileup = 0.25;
35  // Next are vaules I used when I first sent data to Lanny.
36  // double fPileup[] = {0.0, 0.0, 0.2, 0.3, 0.3, 0.2, 0.1, 0.05, 0.0, 0.0, 0.0};
37  } else if (!strcmp(dataset,"AuAu62GeV")) {
38  int nCent = 11;
39  TFile *data1 = new TFile("AuAu62GeV_11c_pair.root");
40  TFile *data2 = new TFile("AuAu62GeV_11c_pair_noPileupCuts.root");
41  const char* centName[] = {"90-100", "80-90", "70-80", "60-70", "50-60", "40-50", "30-40", "20-30", "10-20", "5-10", "0-5"};
42  double norm[] = {1.0298, 1.1747, 1.1758, 1.1937, 1.2174, 1.2416, 1.2702, 1.3082, 1.3668, 1.4016, 1.4426};
43  double fPileup = 0.25;
44  } else if (!strcmp(dataset,"CuCu200GeV")) {
45  int nCent = 10;
46  TFile *data1 = new TFile("CuCu200GeV_10c_pair.root");
47  TFile *data2 = new TFile("CuCu200GeV_10c_pair_noPileupCuts.root");
48  const char* centName[] = {"90-100", "80-90", "70-80", "60-70", "50-60", "40-50", "30-40", "20-30", "10-20", "5-10", "0-5"};
49  double fPileup = 0.25;
50  double norm[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
51  } else if (!strcmp(dataset,"CuCu62GeV")) {
52  int nCent = 10;
53  TFile *data1 = new TFile("CuCu62GeV_10c_pair.root");
54  TFile *data2 = new TFile("CuCu62GeV_10c_pair_noPileupCuts.root");
55  const char* centName[] = {"90-100", "80-90", "70-80", "60-70", "50-60", "40-50", "30-40", "20-30", "10-20", "5-10", "0-5"};
56  double fPileup = 0.25;
57  double norm[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
58  }
59 }
60 
61 // The details of getting pointers to histograms depends on the cutBin mode.
62 // Here we assume mode 3.
63 
64 const char* binName[] = {"all","soft","neck","hard","softHard"};
65 const char* chargeName[] = {"LS", "US", "CD", "CI"};
66 const char* chargeType[] = {"_PP_", "_PM_", "_MP_", "_MM_"};
67 
68 TH2D *dedpCut[nCent];
69 TH2D *dedpNoCut[nCent];
70 TH2D *dedpNoPileup[nCent];
71 TH2D *dedpPileup[nCent];
72 TH2D *ptdedp[nCent];
73 
74  // Get pointers to histograms.
75 {
76  for (int ic=0;ic<nCent;ic++) {
77  data1->cd();
78  TString name(binName[ibin]);
79  name += "_NDEtaDPhi_"; name += chargeName[icharge]; name += "_"; name += ic;
80  dedpCut[ic] = (TH2D *) gDirectory->Get(name.Data());
81  dedpCut[ic]->Scale(norm[ic]);
82  TString name(binName[ibin]);
83  name += "_PtDEtaDPhi_"; name += chargeName[icharge]; name += "_"; name += ic;
84  ptdedp[ic] = (TH2D *) gDirectory->Get(name.Data());
85  ptdedp[ic]->Scale(norm[ic]);
86 
87  data2->cd();
88  TString name(binName[ibin]);
89  name += "_NDEtaDPhi_"; name += chargeName[icharge]; name += "_"; name += ic;
90  dedpNoCut[ic] = (TH2D *) gDirectory->Get(name.Data());
91  dedpNoCut[ic]->Scale(norm[ic]);
92 
93  dedpNoPileup[ic] = NULL;
94  dedpPileup[ic] = NULL;
95  }
96 }
97 
98  // Extrapolate to no pileup
99  // Since Cut and noCut are mostly the same events they don't have independent errors.
100  // Root propagates errors to noPileup as if they are independent. We need to reset those.
101 {
102  for (int ic=0;ic<nCent;ic++) {
103  if (dedpNoPileup[ic]) {
104  delete dedpNoPileup[ic];
105  dedpNoPileup[ic] = NULL;
106  }
107  if (dedpPileup[ic]) {
108  delete dedpPileup[ic];
109  dedpPileup[ic] = NULL;
110  }
111  dedpNoPileup[ic] = (TH2D *) dedpCut[ic]->Clone();
112  dedpPileup[ic] = (TH2D *) dedpCut[ic]->Clone();
113  dedpNoPileup[ic]->Reset();
114  dedpPileup[ic]->Reset();
115 
116  dedpNoPileup[ic]->Add(dedpCut[ic],dedpNoCut[ic],1,-fPileup);
117  dedpPileup[ic]->Add(dedpNoCut[ic],dedpCut[ic],1,-1);
118  dedpNoPileup[ic]->Scale(1/(1-fPileup));
119  dedpPileup[ic]->Scale(1/(1-fPileup));
120  for (int ix=1;ix<=dedpNoPileup[ic]->GetNbinsX();ix++) {
121  for (int iy=1;iy<=dedpNoPileup[ic]->GetNbinsY();iy++) {
122  double err = dedpCut[ic]->GetBinError(ix,iy);
123  dedpNoPileup[ic]->SetBinError(ix,iy,err);
124  dedpPileup[ic]->SetBinError(ix,iy,err);
125  }
126  }
127  }
128 }
129 
130 // For CuCu the bin at (0,0) is screwy, LS histogram is just missing here.
131 // Set error to large number (or could set it to 0 to remove from fit.
132 {
133  if (!strcmp(dataset,"CuCu200GeV") || !strcmp(dataset,"CuCu62GeV")) {
134  for (int ic=0;ic<nCent;ic++) {
135  dedpNoPileup[ic]->SetBinError(13,7,1);
136  dedpPileup[ic]->SetBinError(13,7,1);
137  }
138  }
139 }
140 
141 //-------------------------------------------------
142 // Data has been read.
143 // Create fit functions.
144 .L fitCos.C++
145 
146 TAxis *x = dedpCut[3][0]->GetXaxis();
147 TAxis *y = dedpCut[3][0]->GetYaxis();
148 int nx = x->GetNbins();
149 double xmin = x->GetXmin();
150 double xmax = x->GetXmax();
151 int ny = y->GetNbins();
152 double ymin = y->GetXmin();
153 double ymax = y->GetXmax();
154 double pi = 3.1415926;
155 
156 TF2 *fitFunc[nCent];
157 {
158  char buffer[1024];
159  for (int ic=0;ic<nCent;ic++) {
160  sprintf(buffer,"fit_bin%i_charge%i_Cent%i",ibin,icharge,ic);
161  fitFunc[ic] = new TF2(buffer,fitCos,xmin,xmax,ymin,ymax,nParams);
162  }
163 }
164 
165 // For covariance matrix.
166 double matrix[nCent][nParams][nParams];
167 {
168  for (int ic=0;ic<nCent;ic++) {
169  for (int ip=0;ip<nParams;ip++) {
170  for (int jp=0;jp<nParams;jp++) {
171  matrix[ic][ip][jp] = 0;
172  }
173  }
174  }
175 }
176 
177 //-------------------------------------------------
178 // Read parameters from previous fit. Set lastFits to NULL if we don't have previous
179 // fit (or the format for the previous fit is different.)
180 //char *lastFits = NULL;
181 
182 char lastFits[1024];
183 sprintf(lastFits,"fitResults/cosEta_fitParams_%s_%ic.txt",dataset,nCent);
184 .L readParams.C++
185 {
186  if (lastFits) {
187  float lastParams[nCent][nParams], lastErrors[nCent][nParams];
188  readParams(lastFits,&lastParams[0][0],&lastErrors[0][0],nCent,nParams);
189  }
190 }
191 
192 // Set initial values and limits on parameters.
193 {
194  char *parName[] = {"offset", "cos(phi)", "cos(2phi)", "jetAmp", "jetSigX", "jetSigY", "expAmp", "expX", "expY", "cos(etaD)", "cos(2etaD)", "cos(3etaD)", "cos(4etaD)"};
195  double parVal[] = {0.0, -0.03, 0.03, 0.15,0.6,0.6, 0.1,0.2,0.2, -0.1, 0.05, -0.02, 0.01};
196  for (int ic=0;ic<nCent;ic++) {
197  for (int ip=0;ip<nParams;ip++) {
198  fitFunc[ic]->SetParName(ip,parName[ip]);
199  if (lastFits) {
200  fitFunc[ic]->SetParameter(ip,lastParams[ic][ip]);
201  } else {
202  fitFunc[ic]->SetParameter(ip,parVal[ip]);
203  }
204  }
205  }
206 }
207 {
208  for (int ic=0;ic<nCent;ic++) {
209  fitFunc[ic]->SetParLimits( 0,-2.0,2.0);
210  fitFunc[ic]->SetParLimits( 1,-1.6,1.0);
211  fitFunc[ic]->SetParLimits( 2,-1.0,1.0);
212  fitFunc[ic]->SetParLimits( 3,0.0,2.0);
213  fitFunc[ic]->SetParLimits( 4,0.1,5.0);
214  fitFunc[ic]->SetParLimits( 5,0.1,5.0);
215  fitFunc[ic]->SetParLimits( 6,0.0,25.0);
216  fitFunc[ic]->SetParLimits( 7,1e-9,0.95);
217  fitFunc[ic]->SetParLimits( 8,1e-9,0.95);
218  fitFunc[ic]->SetParLimits( 9,-1.0,1.0);
219  fitFunc[ic]->SetParLimits(10,-1.0,1.0);
220  fitFunc[ic]->SetParLimits(11,-1.0,1.0);
221  fitFunc[ic]->SetParLimits(12,-1.0,1.0);
222  }
223 }
224 
225 // Initial fits. No MINOS errors so fit should be fairly quick.
226 {
227  for (int ic=0;ic<nCent;ic++) {
228  reject = kTRUE;
229  printf(">>Starting fit to centrality %i, bin %i, charge %i\n",ic,ibin,icharge);
230  sprintf(buffer,"fit_bin%i_charge%i_Cent%i",ibin,icharge,ic);
231  dedpNoPileup[ic]->Fit(buffer,"O");
232  }
233 }
234 
235 // Free selected paramters (may make a difference in errors).
236 {
237  for (int ic=0;ic<nCent;ic++) {
238  fitFunc[ic]->ReleaseParameter( 0);
239  fitFunc[ic]->ReleaseParameter( 1);
240  fitFunc[ic]->ReleaseParameter( 2);
241  fitFunc[ic]->ReleaseParameter( 3);
242  fitFunc[ic]->ReleaseParameter( 4);
243  fitFunc[ic]->ReleaseParameter( 5);
244  fitFunc[ic]->ReleaseParameter( 6);
245  fitFunc[ic]->ReleaseParameter( 7);
246  fitFunc[ic]->ReleaseParameter( 8);
247  fitFunc[ic]->ReleaseParameter( 9);
248  fitFunc[ic]->ReleaseParameter(10);
249  }
250 }
251 
252 
253 // Redo fit asking for MINOS errors.
254 {
255  for (int ic=0;ic<nCent;ic++) {
256  reject = kTRUE;
257  printf(">>Starting fit to centrality %i, bin %i, charge %i\n",ic,ibin,icharge);
258  sprintf(buffer,"fit_bin%i_charge%i_Cent%i",ibin,icharge,ic);
259  dedpNoPileup[ic]->Fit(buffer,"EO");
260  gMinuit->mnemat(&matrix[ic][0][0],nParams);
261  }
262 }
263 
264 //-------------------------------------------------
265 // Extract fit parameters, write to file.
266 reject = kFALSE;
267 double param[nCent][nParams];
268 double error[nCent][nParams];
269 double chi2[nCent];
270 {
271  for (int ic=0;ic<nCent;ic++) {
272  fitFunc[ic]->GetParameters(param[ic]);
273  chi2[ic] = fitFunc[ic]->GetChisquare();
274  for (int ip=0;ip<nParams;ip++) {
275  error[ic][ip] = fitFunc[ic]->GetParError(ip);
276  }
277  }
278 }
279 char *paramNames[] = {" chi2 ", " offset", " cos(phi)", " cos(2phi)", " jetAmp", " jetSigX", " jetSigY", " expAmp", " expSigX", " expSigY", " cos(etaD)", " cos(2etaD)", " cos(3etaD)", " cos(4etaD)"};
280 .L writeParams.C++
281 char newFits[1024];
282 sprintf(newFits,"fitResults/cosEta_fitParams_%s_%ic.txt",dataset,nCent);
283 FILE *fPar = fopen(newFits,"w");
284 writeParams(fPar,ibin,icharge,paramNames,chi2,&param[0][0],&error[0][0],nCent,nParams);
285 fclose(fPar);
286 
287 // Write parameters in format appropriate for creating root plots.
288 char *pNames[] = {"offset", "cosPhi", "cos2Phi", "MjAmp", "MjEta", "MjPhi", "expAmp", "expEta", "expPhi", "cos_etaD", "cos_2etaD", "cos_3etaD", "cos_4etaD"};
289 .L writePlotCode.C++
290 char plotCode[1024];
291 sprintf(plotCode,"fitResults/cos_plotCode_%s_%ic.C",dataset,nCent);
292 char tag[1024];
293 sprintf(tag,"cos%s_%i_",dataset,nCent);
294 writePlotCode(plotCode,tag,pNames,chi2,&param[0][0],&error[0][0],nCent,nParams);
295 
296 // Write normalized covariance matrices.
297 // Print covariance matrices.
298 // Seems that Minuit omits fixed parameters.
299 char *cNames[] = {" offset", " cosPhi", " cos2Phi", " mjAmp", " mjEta", " mjPhi", " eAmp", " eEta", " ePhi", " cEta", " c2Eta", " c3Eta", " c4Eta"};
300 .L writeCovariance.C++
301 char covariance[1024];
302 sprintf(covariance,"fitResults/cos_covariance_%s_%ic.txt",dataset,nCent);
303 FILE *fPar = fopen(covariance,"w");
304 writeCovariance(fPar,ibin,icharge,cNames,error[0],matrix[0][0],nCent,nParams);
305 fclose(fPar);
306 
307 
308 
309 //==================================================================================================
310 // Most of rest of this code is for drawing histograms, fits, residuals, etc.
311 // Basically checking quality of fits.
312 
313 TCanvas* c1=new TCanvas("c1");
314 gStyle->SetPalette(1); // set up the colors
315 c1->Clear();
316 //gStyle->SetOptTitle(0);
317 gStyle->SetOptStat(0);
318 c1->SetWindowSize(800,600);
319 c1->Divide(4,3);
320 
321 //-------------------------------------------------
322 // Drawing fitFunction often does not get the sharp exponential.
323 // Copy to histograms.
324 TH2D *fitHist[nCent];
325 reject = kFALSE;
326 {
327  for (int ic=0;ic<nCent;ic++) {
328  fitHist[ic] = (TH2D *) dedpNoPileup[ic]->Clone();
329  fitHist[ic]->Reset();
330  fitHist[ic]->Add(fitFunc[ic]);
331  }
332 }
333 
334 
335 reject = kFALSE;
336 {
337  for (int ic=0;ic<nCent;ic++) {
338  c1->cd(1+ic);
339  gPad->SetPhi(30);
340  gPad->SetTheta(30);
341  fitHist[ic]->Draw("surf1,hist");
342  }
343 }
344 
345 //-------------------------------------------------
346 {
347  for (int ic=0;ic<nCent;ic++) {
348  c1->cd(ic+1);
349  gPad->SetPhi(30);
350  gPad->SetTheta(30);
351  dedpNoPileup[ic]->Draw("surf1,hist");
352  }
353 }
354 
355 //-------------------------------------------------
356 // Subtract fit from data to get residual.
357 TH2D *residual[nCent];
358 {
359  for (int ic=0;ic<nCent;ic++) {
360  residual[ic] = (TH2D *) dedpNoPileup[ic]->Clone();
361  residual[ic]->Add(fitHist[ic],-1);
362  double max = dedpNoPileup[ic]->GetMaximum();
363  double min = dedpNoPileup[ic]->GetMinimum();
364  residual[ic]->SetMaximum(0.05*(max-min));
365  residual[ic]->SetMinimum(-0.05*(max-min));
366  }
367 }
368 
369 {
370  for (int ic=0;ic<nCent;ic++) {
371  c1->cd(ic+1);
372  gPad->SetPhi(30);
373  gPad->SetTheta(50);
374  residual[ic]->Draw("surf1,hist");
375  }
376 }
377 
378 //-------------------------------------------------
379 // Divide residual by error to get chi squared map.
380 TH2D *chiSqMap[nCent];
381 reject = kFALSE;
382 {
383  for (int ic=0;ic<nCent;ic++) {
384  chiSqMap[ic] = NULL;
385  }
386 }
387 {
388  reject = kFALSE;
389  for (int ic=0;ic<nCent;ic++) {
390  if (chiSqMap[ic]) {
391  delete chiSqMap[ic];
392  chiSqMap[ic] = NULL;
393  }
394  chiSqMap[ic] = (TH2D *) dedpNoPileup[ic]->Clone();
395  chiSqMap[ic]->Add(fitFunc[ic],-1);
396  for (int ix=1;ix<=chiSqMap[ic]->GetNbinsX();ix++) {
397  for (int iy=1;iy<=chiSqMap[ic]->GetNbinsY();iy++) {
398  double den = dedpNoPileup[ic]->GetBinError(ix,iy);
399  if (den != 0) {
400  double res = chiSqMap[ic]->GetBinContent(ix,iy)/den;
401  chiSqMap[ic]->SetBinContent(ix,iy,res);
402  }
403  }
404  }
405  }
406 }
407 
408 {
409  for (int ic=0;ic<nCent;ic++) {
410  c1->cd(1+ic);
411  gPad->SetPhi(30);
412  gPad->SetTheta(50);
413 // chiSqMap[ic][ibin]->SetMaximum(max[icharge]);
414  chiSqMap[ic]->Draw("colz,hist");
415  }
416 }
417 
418 //-------------------------------------------------
419 // Plot errors
420 TH2D *hErr[nCent];
421 {
422  for (int ic=0;ic<nCent;ic++) {
423  hErr[ic] = (TH2D *) dedpNoPileup[ic]->Clone();
424  for (int ix=1;ix<=hErr[ic]->GetNbinsX();ix++) {
425  for (int iy=1;iy<=hErr[ic]->GetNbinsY();iy++) {
426  double e = dedpNoPileup[ic]->GetBinError(ix,iy);
427  hErr[ic]->SetBinContent(ix,iy,e);
428  }
429  }
430  }
431 }
432 
433 {
434  for (int ic=0;ic<nCent;ic++) {
435  c1->cd(ic+1);
436  gPad->SetPhi(30);
437  gPad->SetTheta(30);
438  hErr[ic]->Draw("surf1,hist");
439  }
440 }
441 
442 
443 
444 //-------------------------------------------------
445 // Write data and errors as text files.
446 {
447  char buffer[1024];
448  for (int ic=0;ic<nCent;ic++) {
449  sprintf(buffer,"bin%i_CIDEtaDPhi_softHard.txt",ic);
450  FILE *fAmp = fopen(buffer,"w");
451  sprintf(buffer,"bin%i_CIDEtaDPhi_Errors_softHard.txt",ic);
452  FILE *fErr = fopen(buffer,"w");
453  double eta, phi, a, e;
454  for (int ix=1;ix<=dedpNoPileup[ic]->GetNbinsX();ix++) {
455  for (int iy=1;iy<=dedpNoPileup[ic]->GetNbinsY();iy++) {
456  eta = dedpNoPileup[ic]->GetXaxis()->GetBinCenter(ix);
457  phi = dedpNoPileup[ic]->GetYaxis()->GetBinCenter(iy);
458  a = dedpNoPileup[ic]->GetBinContent(ix,iy);
459  e = dedpNoPileup[ic]->GetBinError(ix,iy);
460  fprintf(fAmp," %7.4f %7.4f %7.4f\n",eta,phi,a);
461  fprintf(fErr," %7.4f %7.4f %9.6f\n",eta,phi,e);
462  }
463  }
464  fclose(fAmp);
465  fclose(fErr);
466  }
467 }
468 
469 //-------------------------------------------------
470 // Sometimes we want to look at particular components.
471 TF2 *fitComponent[nCent][9];
472 char buffer[1024];
473 {
474  for (int ic=0;ic<nCent;ic++) {
475  sprintf(buffer,"fit_Cent%i_Component0",ic);
476  fitComponent[ic][0] = new TF2(buffer,"[0]",xmin,xmax,ymin,ymax);
477  fitComponent[ic][0]->SetNpx(nx);
478  fitComponent[ic][0]->SetNpy(ny);
479  sprintf(buffer,"fit_Cent%i_Component1",ic);
480  fitComponent[ic][1] = new TF2(buffer,"[0]*cos(y)",xmin,xmax,ymin,ymax);
481  fitComponent[ic][1]->SetNpx(nx);
482  fitComponent[ic][1]->SetNpy(ny);
483  sprintf(buffer,"fit_Cent%i_Component2",ic);
484  fitComponent[ic][2] = new TF2(buffer,"[0]*cos(2*y)",xmin,xmax,ymin,ymax);
485  fitComponent[ic][2]->SetNpx(nx);
486  fitComponent[ic][2]->SetNpy(ny);
487  sprintf(buffer,"fit_Cent%i_Component3",ic);
488  fitComponent[ic][3] = new TF2(buffer,"[0]*exp(-0.5*((x/[1])^2+(y/[2])^2))+[0]*exp(-0.5*((x/[1])^2+((y-6.283)/[2])^2))",xmin,xmax,ymin,ymax);
489  fitComponent[ic][3]->SetNpx(nx);
490  fitComponent[ic][3]->SetNpy(ny);
491  sprintf(buffer,"fit_Cent%i_Component4",ic);
492  fitComponent[ic][4] = new TF2(buffer,"[0]*exp(-sqrt((x/[1])^2+(y/[2])^2))",xmin,xmax,ymin,ymax);
493  fitComponent[ic][4]->SetNpx(nx);
494  fitComponent[ic][4]->SetNpy(ny);
495  sprintf(buffer,"fit_Cent%i_Component5",ic);
496  fitComponent[ic][5] = new TF2(buffer,"[0]*cos(x)",xmin,xmax,ymin,ymax);
497  fitComponent[ic][5]->SetNpx(nx);
498  fitComponent[ic][5]->SetNpy(ny);
499  sprintf(buffer,"fit_Cent%i_Component6",ic);
500  fitComponent[ic][6] = new TF2(buffer,"[0]*cos(2*x)",xmin,xmax,ymin,ymax);
501  fitComponent[ic][6]->SetNpx(nx);
502  fitComponent[ic][6]->SetNpy(ny);
503  sprintf(buffer,"fit_Cent%i_Component7",ic);
504  fitComponent[ic][7] = new TF2(buffer,"[0]*cos(3*x)",xmin,xmax,ymin,ymax);
505  fitComponent[ic][7]->SetNpx(nx);
506  fitComponent[ic][7]->SetNpy(ny);
507  sprintf(buffer,"fit_Cent%i_Component8",ic);
508  fitComponent[ic][8] = new TF2(buffer,"[0]*cos(4*x)",xmin,xmax,ymin,ymax);
509  fitComponent[ic][8]->SetNpx(nx);
510  fitComponent[ic][8]->SetNpy(ny);
511  }
512 }
513 {
514  char buffer[1024];
515  for (int ic=0;ic<nCent;ic++) {
516  fitComponent[ic][0]->SetParameter(0,param[ic][0]);
517 
518  fitComponent[ic][1]->SetParameter(0,param[ic][1]);
519 
520  fitComponent[ic][2]->SetParameter(0,param[ic][2]);
521 
522  fitComponent[ic][3]->SetParameter(0,param[ic][3]);
523  fitComponent[ic][3]->SetParameter(1,param[ic][4]);
524  fitComponent[ic][3]->SetParameter(2,param[ic][5]);
525 
526  fitComponent[ic][4]->SetParameter(0,param[ic][6]);
527  fitComponent[ic][4]->SetParameter(1,param[ic][7]);
528  fitComponent[ic][4]->SetParameter(2,param[ic][8]);
529 
530  fitComponent[ic][5]->SetParameter(0,param[ic][9]);
531  fitComponent[ic][6]->SetParameter(0,param[ic][10]);
532  fitComponent[ic][7]->SetParameter(0,param[ic][11]);
533  fitComponent[ic][8]->SetParameter(0,param[ic][12]);
534  }
535 }
536 {
537  for (int ic=0;ic<nCent;ic++) {
538  c1->cd(ic+1);
539  gPad->SetPhi(30);
540  gPad->SetTheta(30);
541  fitComponent[ic][4]->Draw("surf1,hist");
542  }
543 }
544 
545 TH2D *hSub[nCent];
546 {
547  for (int ic=0;ic<nCent;ic++) {
548  hSub[ic] = (TH2D *) dedpNoPileup[ic]->Clone();
549  hSub[ic]->Reset();
550  hSub[ic]->Add(fitComponent[ic][5],1);
551  hSub[ic]->Add(fitComponent[ic][6],1);
552  hSub[ic]->Add(fitComponent[ic][7],1);
553  hSub[ic]->Add(fitComponent[ic][8],1);
554  }
555 }
556 {
557  for (int ic=0;ic<nCent;ic++) {
558  hSub[ic]->Add(fitComponent[ic][3],-1);
559  }
560 }
561 {
562  for (int ic=0;ic<nCent;ic++) {
563  c1->cd(ic+1);
564  gPad->SetPhi(30);
565  gPad->SetTheta(30);
566  hSub[ic]->Draw("surf1,hist");
567  }
568 }
569 
570 
571 //-------------------------------------------------
572 gStyle->SetTitleBorderSize(0);
573 gPad->SetLeftMargin(0.16);
574 gROOT->SetStyle("Plain");
575 gStyle->SetOptStat(0);
576 
577 hSub[6]->SetTitle("Same-side peak, 28-38% centrality")
578 TAxis *x = hSub[6]->GetXaxis();
579 x->SetTitleSize(0.07);
580 x->SetTitleOffset(1.0);
581 x->SetNdivisions(505);
582 x->SetLabelSize(0.05);
583 x->SetTitle("#eta_{#Delta}");
584 TAxis *y = hSub[6]->GetYaxis();
585 y->SetTitleSize(0.07);
586 y->SetTitleOffset(1.0);
587 y->SetNdivisions(505);
588 y->SetLabelSize(0.05);
589 y->SetTitle("#phi_{#Delta}");
590 TAxis *z = hSub[6]->GetZaxis();
591 z->SetTitleSize(0.07);
592 z->SetTitleOffset(1.0);
593 z->SetNdivisions(505);
594 z->SetLabelSize(0.05);
595 z->SetTitle("#Delta#rho/#sqrt{#rho_{ref}}");
596 
597 hSub[4]->SetTitle("Same-side peak, 46-55% centrality")
598 TAxis *x = hSub[4]->GetXaxis();
599 x->SetTitleSize(0.07);
600 x->SetTitleOffset(1.0);
601 x->SetNdivisions(505);
602 x->SetLabelSize(0.05);
603 x->SetTitle("#eta_{#Delta}");
604 TAxis *y = hSub[4]->GetYaxis();
605 y->SetTitleSize(0.07);
606 y->SetTitleOffset(1.0);
607 y->SetNdivisions(505);
608 y->SetLabelSize(0.05);
609 y->SetTitle("#phi_{#Delta}");
610 TAxis *z = hSub[4]->GetZaxis();
611 z->SetTitleSize(0.07);
612 z->SetTitleOffset(1.0);
613 z->SetNdivisions(505);
614 z->SetLabelSize(0.05);
615 z->SetTitle("#Delta#rho/#sqrt{#rho_{ref}}");
Generated by   doxygen 1.8.5