added LHS gaus fit (Marta)
[u/mrichter/AliRoot.git] / PWG4 / macros / jetspectrum / PlotNote.C
CommitLineData
e2081c22 1#include "TList.h"
2#include "TColor.h"
3#include "TArrow.h"
4#include "TROOT.h"
5#include "TLatex.h"
6#include "TSystem.h"
7#include "TGraph.h"
8#include "TGraphErrors.h"
9#include "TFile.h"
10#include "TLegend.h"
11#include "TStyle.h"
12#include "TCanvas.h"
13#include "THnSparse.h"
14#include "TDirectoryFile.h"
15#include "TH1F.h"
16#include "TH2F.h"
17#include "TH3F.h"
18#include "TMath.h"
19#include "TF1.h"
20#include "TASImage.h"
21#include "TPaveText.h"
22#include "TParameter.h"
23#include "TPaveStats.h"
24
25
26#include "AliAnalysisHelperJetTasks.h"
27#include "AliAnalysisTaskJetServices.h"
28#include "AliAnalysisTaskJetSpectrum2.h"
29#include <fstream>
30#include <iostream>
31
32#include "Normalize2d.C"
33#include "ALICEWorkInProgress.C"
34
35using namespace std;
36
37Int_t PlotSpectraPbPb();
38Int_t PlotJetBFluctuations();
39Int_t PlotJetQA();
40// helpers
41
42
43enum {kMaxJets = 4};
44void HistsFromSingleJets(char* cMask,TList *list,Int_t nJets,TH1F** h,Int_t iRebin = 1,int iFirst = 0);
45void HistsFromSingleJets(char* cMask,TList *list,Int_t nJets,TH2F** h,Int_t iRebin = 1,int iFirst = 0);
46void ScaleNevent(TH1 *h,TFile *fIn,Int_t iCond = 5,Int_t iMC = 0,Int_t iCen = 0);
47void ScaleH1(TH1* h,char* cX = "",char* cY = "",Float_t fScale = 1,Bool_t bWidth = true);
48void SetStyleH1(TH1 *h,Int_t color = 0,Int_t fillStyle = 0,Int_t markerStyle = 0);
49void ScaleH2(TH2* h,char* cX = "",char* cY = "",char* cZ = "",Float_t fScale = 1,Bool_t bWidth = true);
50TH1F* GetPythia8Spectrum(const char *cHist,const char *cPythia,Float_t deta,Int_t iRebin = 1,Int_t iMask = 0);
51const char *cPrintMask = "110116_%s.eps";
52void set_plot_style();
53TCanvas *cTmp = 0;
54
55TH1F* GetPythia8Spectrum(const char *cHist,const char *cPythia,Float_t deta,Int_t iRebin,Int_t iMask){
56 TDirectory *opwd = gDirectory;
57 TFile *fPythia = (TFile*)gROOT->GetListOfFiles()->FindObject(cPythia);
58
59 if(!fPythia) fPythia = TFile::Open(cPythia);
60 opwd->cd();
61
62 static int iCount = 0;
63 TH1F* hist = 0;
64 if(iMask<=0){
65 hist = (TH1F*)fPythia->Get(cHist);
66 if(hist)hist = (TH1F*)hist->Clone(Form("%s_%d",hist->GetName(),iCount));
67 }
68 else{
69 for(int i = 0;i < iMask;i++){
70 TH1F *hTmp = (TH1F*)fPythia->Get(Form(cHist,i));
71 if(!hTmp){
72 Printf("%s not found",Form(cHist,i));
73 return 0;
74 }
75 if(!hist)hist = (TH1F*)hTmp->Clone(Form("%s_%d_%d",hTmp->GetName(),iMask,iCount));
76 else hist->Add(hTmp);
77 }
78 }
79
80 if(!hist){
81 Printf("%s not found",cHist);
82 return 0;
83 }
84 // fetch the cross section
85 TH1F *hxsec = (TH1F*)fPythia->Get("h1Xsec");
86
87 Float_t xsec = hxsec->GetBinContent(1);
88 Printf("%d xsec = %1.3E",__LINE__,xsec);
89 // if(xsec==0)xsec = 40.79; // tmp hack
90 hist->Rebin(iRebin);
91 hist->Scale(1./hist->GetBinWidth(1));
92 // hist->Scale(1./xsec);
93 hist->Scale(1./deta);
94
95 // Double_t xtotal = 71.39; // xtotal for 7 TeV!!
96 Double_t xtotal = 62.01; // xtotal for 7 TeV!!
97
98 // scale with fraction of total NSD cross section
99 hist->Scale(1/xtotal);
100 Printf("%d fraction = %1.3E of total cross section %1.3E",__LINE__,xsec/xtotal,xtotal);
101 hist->SetXTitle("p_{T} (GeV/c)");
102 hist->SetYTitle("1/N dN/dp_{T}dy");
103 hist->SetDirectory(opwd);
104 iCount++;
105 return hist;
106
107}
108
109void PlotNote(){
110
111 PlotJetBFluctuations();
112}
113
114TList *GetList(const char *cFile,const char *cName){
115 TDirectory *opwd = gDirectory;
116 TFile *fIn = (TFile*)gROOT->GetListOfFiles()->FindObject(cFile);
117 TList *list = 0;
118 if(!fIn) fIn = TFile::Open(cFile);
119 if(!fIn)return list;
120 opwd->cd();
121
122 TDirectory *dir = (TDirectory*)fIn->Get(Form("PWG4_%s",cName));
123 if(!dir){
124 Printf("GetList: Directory PWG4_%s not found",cName);
125 return list;
126 }
127 list = (TList*)dir->Get(Form("pwg4%s",cName));
128 if(!list){
129 Printf("GetList: list pwg4%s not found",cName);
130 return list;
131 }
132 return list;
133
134}
135
136Int_t PlotJetQA(){
137
138 set_plot_style();
139 gStyle->SetTitleYSize(0.7*gStyle->GetTitleYSize());
140 gStyle->SetTitleXSize(0.7*gStyle->GetTitleXSize());
141 gStyle->SetTitleXOffset(1.2*gStyle->GetTitleXOffset());
142 gStyle->SetTitleYOffset(1.5*gStyle->GetTitleYOffset());
143 gStyle->SetPadRightMargin(1.1*gStyle->GetPadRightMargin());
144 gStyle->SetPadBottomMargin(0.7*gStyle->GetPadBottomMargin());
145
146
147 Bool_t isPP = false;
148
149 TString printType = "eps";
150 const Int_t nCen = 5;
151 TString sCen[nCen] = {" 0-80"," 0-10%%","10-30","30-50","50-80"};
152 const Int_t nJets = 3;
153
154
155 TString cAdd = "Rec"; // decide whcih type to take Rec RecFull Gen GenFull
156 TString cAddFull = "RecFull"; // decide whcih type to take Rec RecFull Gen GenFull
157
158 TCanvas *c1 = new TCanvas();
159
160 TString sFile = "~/alice/data/analysis/train_maf/output_110216.root";
161 TString sinputName = "spec2_clustersAOD_ANTIKT04_B2_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000";
162 TString sinputJet = "%% Pb+Pb Anti-k_{T} R = 0.4 (B2)";
163 TString sinputPrint = "PbPb_antikt04_B2";
164 sinputName = "spec2_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B0_Filter00256_Cut00150_Skip00_0000000000";cAdd = "Gen";cAddFull = "GenFull";sinputJet = "%% Pb+Pb Anti-k_{T} R = 0.4 (B0)";sinputPrint = "PbPb_antikt04_B0";
165
166 if(isPP){
167 sinputName = "spec2_clustersAOD_ANTIKT04_B0_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B0_Filter00256_Cut00150_Skip00_0000000000";sinputJet = "p+p #sqrt{s} = 7 TeV Anti-k_{T} R = 0.4 (B0)";sinputPrint = "pp_antikt04_B0";
168 sCen[0] = "";
169 sFile = "~/alice/jets/macros/batch/output.root";
170 }
171
172
173 TFile *fIn = TFile::Open(sFile.Data());
174
175
176 TDirectory *dir = (TDirectory*)fIn->Get("PWG4_services");
177 TList *list1 = (TList*)dir->Get("pwg4serv");
178
179 TH2F* hTriggerCount =(TH2F*)list1->FindObject("fh2TriggerCount");
180
181 for(int ic = 0;ic < (isPP?1:nCen);ic++){
182
183 TString cName = Form("PWG4_%s_Class%02d",sinputName.Data(),ic);
184 TDirectory *inputDir = (TDirectory*)fIn->Get(cName.Data());
185
186 Float_t fNevents = hTriggerCount->GetBinContent(hTriggerCount->FindBin(ic,AliAnalysisTaskJetServices::kSelected));
187
188 Printf(" %s: %10d events",sCen[ic].Data(),(Int_t)fNevents);
189
190 if(!inputDir){
191 Printf("Dir %s not found",cName.Data());
192 continue;
193 }
194 cName = Form("pwg4%s_Class%02d",sinputName.Data(),ic);
195 TList *inputList = (TList*)inputDir->Get(cName.Data());
196 if(!inputList){
197 Printf("List %s not found",cName.Data());
198 continue;
199 }
200
201 TH2F *h2EtaPtLead[nJets+2] = {0,};
202 TH2F *h2EtaPtLeadFull[nJets+2] = {0,};
203 TH2F *h2EtaPtAll[nJets+2] = {0,};
204 HistsFromSingleJets(Form("fh2EtaPt%s_j%%d",cAdd.Data()),inputList,1,h2EtaPtLead);
205 HistsFromSingleJets(Form("fh2EtaPt%s_j%%d",cAddFull.Data()),inputList,1,h2EtaPtLeadFull);
206 HistsFromSingleJets(Form("fh2EtaPt%s_j%%d",cAddFull.Data()),inputList,nJets+1,h2EtaPtAll,1,nJets);
207
208 c1->SetLogz();
209 for(int i = 0; i< nJets+2;i++)Printf("Lead > %d %p",i,h2EtaPtLead[i]);
210 for(int i = 0; i< nJets+2;i++)Printf("LeadFull > %d %p",i,h2EtaPtLeadFull[i]);
211 for(int i = 0; i< nJets+2;i++)Printf("All > %d %p",i,h2EtaPtAll[i]);
212 h2EtaPtLead[0]->Draw("colz");
213 c1->Modified();
214 c1->Update();
215 c1->SaveAs(Form("%s_%s_Cen%02d_EtaPtLead2D.%s",sinputPrint.Data(),cAdd.Data(),ic,printType.Data()));
216 if(!gROOT->IsBatch())
217 if(getchar()=='q')return 1;
218
219 Float_t fStep = 20;
220 Float_t fStepLo = 0;
221 Float_t fStepUp = fStepLo+fStep;
222 Float_t fMax = 160;
223
224 while(fStepUp<fMax){
225 Int_t iBinLo = h2EtaPtLead[0]->GetYaxis()->FindBin(fStepLo);
226 Int_t iBinUp = h2EtaPtLead[0]->GetYaxis()->FindBin(fStepUp)-1;
227 TH1D *hProjLead = h2EtaPtLead[0]->ProjectionX(Form("hEtaPtLead_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
228 TH1D *hProjLeadFull = h2EtaPtLeadFull[0]->ProjectionX(Form("hEtaPtLeadFull_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
229 TH1D *hProjAll = h2EtaPtAll[nJets+1]->ProjectionX(Form("hEtaPtAll_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
230
231 ScaleH1(hProjLead,"#eta","1/N_{evt} dN_{jet}/d#eta",1./fNevents,true);
232 SetStyleH1(hProjLead,kBlue,0,kFullCircle);
233 ScaleH1(hProjLeadFull,"#eta","1/N_{evt} dN_{jet}/d#eta",1./fNevents,true);
234 SetStyleH1(hProjLeadFull,kRed,0,kFullCircle);
235 ScaleH1(hProjAll,"#eta","1/N_{evt} dN_{jet}/d#eta",1./fNevents,true);
236 SetStyleH1(hProjAll,kGray,1001,kFullCircle);
237
238
239 hProjAll->DrawCopy("hist");
240 hProjLead->DrawCopy("histsame");
241 hProjLeadFull->DrawCopy("histsame");
242
243 Float_t ptUp = fStepUp-h2EtaPtLead[0]->GetYaxis()->GetBinWidth(iBinUp);
244
245 Printf("%3.0f-%3.0f",fStepLo,ptUp);
246 TLatex *txt = 0;
247 txt = new TLatex(-1.1,1.1*hProjAll->GetMaximum(),Form("%s %s p_{T} = %3.0f-%3.0f GeV",sCen[ic].Data(),sinputJet.Data(),fStepLo,ptUp));
248 txt->SetTextSize(gStyle->GetTextSize());
249 txt->Draw();
250 c1->SaveAs(Form("%s_%s_Cen%02d_%03d-%03d_eta.%s",sinputPrint.Data(),cAdd.Data(),ic,(Int_t)fStepLo,(Int_t)ptUp,printType.Data()));
251 c1->Modified();
252 c1->Update();
253 if(!gROOT->IsBatch())
254 if(getchar()=='q')return 1;
255 fStepLo += fStep;
256 fStepUp += fStep;
257 delete hProjAll;
258 delete hProjLead;
259 delete hProjLeadFull;
260 }
261
262 // Phi vs. p_T
263
264
265 TH2F *h2PhiPtLead[nJets+2] = {0,};
266 TH2F *h2PhiPtLeadFull[nJets+2] = {0,};
267 TH2F *h2PhiPtAll[nJets+2] = {0,};
268 HistsFromSingleJets(Form("fh2PhiPt%s_j%%d",cAdd.Data()),inputList,1,h2PhiPtLead);
269 HistsFromSingleJets(Form("fh2PhiPt%s_j%%d",cAddFull.Data()),inputList,1,h2PhiPtLeadFull);
270 HistsFromSingleJets(Form("fh2PhiPt%s_j%%d",cAddFull.Data()),inputList,nJets+1,h2PhiPtAll,1,nJets);
271
272 c1->SetLogz();
273 h2PhiPtLead[0]->Draw("colz");
274 c1->SaveAs(Form("%s_%s_Cen%02d_phiPtLead2D.%s",sinputPrint.Data(),cAdd.Data(),ic,printType.Data()));
275 c1->Modified();
276 c1->Update();
277
278 if(!gROOT->IsBatch())
279 if(getchar()=='q')return 1;
280
281 fStepLo = 0;
282 fStepUp = fStepLo+fStep;
283
284 while(fStepUp<fMax){
285 Int_t iBinLo = h2PhiPtLead[0]->GetYaxis()->FindBin(fStepLo);
286 Int_t iBinUp = h2PhiPtLead[0]->GetYaxis()->FindBin(fStepUp)-1;
287 TH1D *hProjLead = h2PhiPtLead[0]->ProjectionX(Form("hPhiPtLead_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
288 TH1D *hProjLeadFull = h2PhiPtLeadFull[0]->ProjectionX(Form("hPhiPtLeadFull_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
289 TH1D *hProjAll = h2PhiPtAll[nJets+1]->ProjectionX(Form("hPhiPtAll_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
290
291 ScaleH1(hProjLead,"#phi","1/N_{evt} dN_{jet}/d#phi",1./fNevents,true);
292 SetStyleH1(hProjLead,kBlue,0,kFullCircle);
293 ScaleH1(hProjLeadFull,"#phi","1/N_{evt} dN_{jet}/d#phi",1./fNevents,true);
294 SetStyleH1(hProjLeadFull,kRed,0,kFullCircle);
295 ScaleH1(hProjAll,"#phi","1/N_{evt} dN_{jet}/d#phi",1./fNevents,true);
296 SetStyleH1(hProjAll,kGray,1001,kFullCircle);
297
298 hProjAll->SetMinimum(0);
299
300 hProjAll->DrawCopy("hist");
301 hProjLead->DrawCopy("histsame");
302 hProjLeadFull->DrawCopy("histsame");
303
304 Float_t ptUp = fStepUp-h2PhiPtLead[0]->GetYaxis()->GetBinWidth(iBinUp);
305
306 Printf("%3.0f-%3.0f",fStepLo,ptUp);
307 TLatex *txt = 0;
308 txt = new TLatex(0.0,1.1*hProjAll->GetMaximum(),Form("%s %s p_{T} = %3.0f-%3.0f GeV",sCen[ic].Data(),sinputJet.Data(),fStepLo,ptUp));
309 txt->SetTextSize(gStyle->GetTextSize());
310 txt->Draw();
311 c1->SaveAs(Form("%s_%s_Cen%02d_%03d-%03d_phi.%s",sinputPrint.Data(),cAdd.Data(),ic,(Int_t)fStepLo,(Int_t)ptUp,printType.Data()));
312 c1->Modified();
313 c1->Update();
314 if(!gROOT->IsBatch())
315 if(getchar()=='q')return 1;
316 fStepLo += fStep;
317 fStepUp += fStep;
318 delete hProjAll;
319 delete hProjLead;
320 delete hProjLeadFull;
321 }
322
323 // area vs. p_T
324
325
326 TH2F *h2AreaPtLead[nJets+2] = {0,};
327 TH2F *h2AreaPtLeadFull[nJets+2] = {0,};
328 TH2F *h2AreaPtAll[nJets+2] = {0,};
329
330 HistsFromSingleJets(Form("fh2AreaPt%s_j%%d",cAdd.Data()),inputList,1,h2AreaPtLead);
331 HistsFromSingleJets(Form("fh2AreaPt%s_j%%d",cAddFull.Data()),inputList,1,h2AreaPtLeadFull);
332 HistsFromSingleJets(Form("fh2AreaPt%s_j%%d",cAddFull.Data()),inputList,nJets+1,h2AreaPtAll,1,nJets);
333
334 c1->SetLogz();
335 h2AreaPtLead[0]->Draw("colz");
336 c1->Modified();
337 c1->Update();
338 c1->SaveAs(Form("%s_%s_Cen%02d_areaPtLead2D.%s",sinputPrint.Data(),cAdd.Data(),ic,printType.Data()));
339 if(!gROOT->IsBatch())
340 if(getchar()=='q')return 1;
341
342 fStepLo = 0;
343 fStepUp = fStepLo+fStep;
344
345 while(fStepUp<fMax){
346 Int_t iBinLo = h2AreaPtLead[0]->GetYaxis()->FindBin(fStepLo);
347 Int_t iBinUp = h2AreaPtLead[0]->GetYaxis()->FindBin(fStepUp)-1;
348 TH1D *hProjLead = h2AreaPtLead[0]->ProjectionX(Form("hAreaPtLead_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
349 TH1D *hProjLeadFull = h2AreaPtLeadFull[0]->ProjectionX(Form("hAreaPtLeadFull_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
350 TH1D *hProjAll = h2AreaPtAll[nJets+1]->ProjectionX(Form("hAreaPtAll_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
351
352 ScaleH1(hProjLead,"Area A","1/N_{evt} dN_{jet}/dA",1./fNevents,true);
353 SetStyleH1(hProjLead,kBlue,0,kFullCircle);
354 ScaleH1(hProjLeadFull,"Area A","1/N_{evt} dN_{jet}/dA",1./fNevents,true);
355 SetStyleH1(hProjLeadFull,kRed,0,kFullCircle);
356 ScaleH1(hProjAll,"Area A","1/N_{evt} dN_{jet}/dA",1./fNevents,true);
357 SetStyleH1(hProjAll,kGray,1001,kFullCircle);
358
359 hProjAll->SetMinimum(0);
360 hProjAll->DrawCopy("hist");
361 hProjLead->DrawCopy("histsame");
362 hProjLeadFull->DrawCopy("histsame");
363
364 Float_t ptUp = fStepUp-h2AreaPtLead[0]->GetYaxis()->GetBinWidth(iBinUp);
365
366 Printf("%3.0f-%3.0f",fStepLo,ptUp);
367 TLatex *txt = 0;
368 txt = new TLatex(0.,1.1*hProjAll->GetMaximum(),Form("%s %s p_{T} = %3.0f-%3.0f GeV",sCen[ic].Data(),sinputJet.Data(),fStepLo,ptUp));
369 txt->SetTextSize(gStyle->GetTextSize());
370 txt->Draw();
371 c1->SaveAs(Form("%s_%s_Cen%02d_%03d-%03d_area.%s",sinputPrint.Data(),cAdd.Data(),ic,(Int_t)fStepLo,(Int_t)ptUp,printType.Data()));
372 c1->Modified();
373 c1->Update();
374 if(!gROOT->IsBatch())
375 if(getchar()=='q')return 1;
376 fStepLo += fStep;
377 fStepUp += fStep;
378 delete hProjAll;
379 delete hProjLead;
380 delete hProjLeadFull;
381 }
382
383 // Area vs eta:
384 TH2F *h2AreaEta[nJets+1] = {0,};
385 HistsFromSingleJets(Form("fh2EtaArea%s_j%%d",cAdd.Data()),inputList,nJets,h2AreaEta);
386 ScaleH2(h2AreaEta[nJets],"area A","#eta","1/N_{evt} dN_{jets}/dAd#eta",1./fNevents,true);
387 c1->SetLogz();
388 h2AreaEta[nJets]->Draw("colz");
389 c1->Modified();
390 c1->Update();
391 c1->SaveAs(Form("%s_%s_Cen%02d_areaEta2D.%s",sinputPrint.Data(),cAdd.Data(),ic,printType.Data()));
392 if(!gROOT->IsBatch())
393 if(getchar()=='q')return 1;
394 }
395
396 return 0;
397}
398
399Int_t PlotSpectraPbPb(){
400
401 // PLot the simple 1D histos from the spectrum task
402 //
403 const Int_t iRebin = 4;
404 const int kMaxFiles = 5;
405
406 Double_t yMin = 0.01;
407 Double_t yMax = 1E+07;
408
409 TString sinputFile[kMaxFiles];
410 TFile* inputFile[kMaxFiles] = {kMaxFiles*0};
411 TString sinputDir[kMaxFiles];
412 TDirectory *inputDir[kMaxFiles] = {kMaxFiles*0};
413 TString sinputList[kMaxFiles];
414 TString sCen[kMaxFiles];
415 TList* inputList[kMaxFiles] = {kMaxFiles*0};
416 // TString sinputLegend[kMaxFiles];
417 TString sinputJet[kMaxFiles];
418 Int_t kColor[kMaxFiles];
419 Float_t nColl[kMaxFiles];
420
421 bool bLogLog = false;
422 // const Int_t kRef = 1; // anti -kT
423
424 //
425 TH1F* hJets[kMaxFiles];
426
427 const int nJets = 2;
428
429 Int_t iJF = 0;
430 sinputFile[iJF] = "~/alice/data/analysis/train/LHC10h_110116/PWG4_JetTasksOutput.root";
431 // sinputFile[iJF] = "~/alice/data/analysis/train_maf/output_110210.root";
432 sinputDir[iJF] = "PWG4_spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class00";
433 sinputList[iJF] = "pwg4spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class00";
434 sinputJet[iJF] = "anti-k_{T} R = 0.4";
435 sCen[iJF] = " 0-80%";
436 kColor[iJF] = kBlack;
437 nColl[iJF] = 453.3;
438 iJF++;
439
440 // sinputFile[iJF] = "../output/PWG4_JetTasksOutput_Merge_051205b_Subtract3_UA1.root";
441 sinputFile[iJF] = "~/alice/data/analysis/train/LHC10h_110116/PWG4_JetTasksOutput.root";
442 // sinputFile[iJF] = "~/alice/data/analysis/train_maf/output_110210.root";
443 sinputDir[iJF] = "PWG4_spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class01";
444 sinputList[iJF] = "pwg4spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class01";
445 sinputJet[iJF] = " 0-10% Anti-k_{T} R = 0.4";
446 sCen[iJF] = " 0-10%";
447 kColor[iJF] = kRed;
448 nColl[iJF] = 1502;
449 iJF++;
450
451 // sinputFile[iJF] = "../output/PWG4_JetTasksOutput_Merge_051205b_Subtract3_UA1.root";
452 sinputFile[iJF] = "~/alice/data/analysis/train/LHC10h_110116/PWG4_JetTasksOutput.root";
453 // sinputFile[iJF] = "~/alice/data/analysis/train_maf/output_110210.root";
454 sinputDir[iJF] = "PWG4_spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class02";
455 sinputList[iJF] = "pwg4spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class02";
456 sinputJet[iJF] = "10-30% Anti-k_{T} R = 0.4";
457 sCen[iJF] = "10-30%";
458 kColor[iJF] = kRed-4;
459 nColl[iJF] = 742.5;
460 iJF++;
461
462 // sinputFile[iJF] = "../output/PWG4_JetTasksOutput_Merge_051205b_Subtract3_UA1.root";
463 sinputFile[iJF] = "~/alice/data/analysis/train/LHC10h_110116/PWG4_JetTasksOutput.root";
464 // sinputFile[iJF] = "~/alice/data/analysis/train_maf/output_110210.root";
465 sinputDir[iJF] = "PWG4_spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class03";
466 sinputList[iJF] = "pwg4spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class03";
467 sinputJet[iJF] = "30-50% Anti-k_{T} R = 0.4";
468 sCen[iJF] = "30-50%";
469 kColor[iJF] = kRed-7;
470 nColl[iJF] = 250;
471 iJF++;
472
473 // sinputFile[iJF] = "../output/PWG4_JetTasksOutput_Merge_051205b_Subtract3_UA1.root";
474 sinputFile[iJF] = "~/alice/data/analysis/train/LHC10h_110116/PWG4_JetTasksOutput.root";
475 // sinputFile[iJF] = "~/alice/data/analysis/train_maf/output_110210.root";
476 sinputDir[iJF] = "PWG4_spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class04";
477 sinputList[iJF] = "pwg4spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class04";
478 sinputJet[iJF] = "Anti-k_{T} R = 0.4";
479 sCen[iJF] = "50-80%";
480 kColor[iJF] = kRed-9;
481 nColl[iJF] = 46.6;
482 iJF++;
483
484
485 for(int iF = 0; iF < kMaxFiles;iF++){
486 if(sinputFile[iF].Length()==0)continue;
487 if(gROOT->GetListOfFiles()->FindObject(sinputFile[iF].Data())){
488 inputFile[iF] = (TFile*)gROOT->GetListOfFiles()->FindObject(sinputFile[iF].Data());
489 }
490 else {
491 inputFile[iF] = TFile::Open(sinputFile[iF].Data());
492 }
493 inputDir[iF] = (TDirectory*)inputFile[iF]->Get(sinputDir[iF].Data());
494 if(!inputDir[iF]){
495 Printf("Dir not found %s",sinputDir[iF].Data());
496 continue;
497 }
498 inputList[iF] = (TList*)inputDir[iF]->Get(sinputList[iF].Data());
499 if(!inputList[iF]){
500 Printf("List not found %s",sinputList[iF].Data());
501 continue;
502 }
503 }
504
505
506
507
508 TCanvas *c1 = new TCanvas("c1","spectra",20,20,800,600);
509 TCanvas *c2 = new TCanvas("c2","ratio",820,20,800,600);
510
511 TLegend *leg1 = new TLegend(0.45,0.65,0.7,0.85);
512 leg1->SetHeader(Form("%s |#eta| < 0.4",sinputJet[0].Data()));
513 leg1->SetFillColor(0);
514 leg1->SetTextFont(gStyle->GetTextFont());
515 leg1->SetBorderSize(0);
516
517 char *cMaskRec = "fh1PtRecIn_j%d";
518
519 c1->SetLogy();
520 if(bLogLog)c1->SetLogx();
521
522 for(int iF = 0; iF < kMaxFiles;iF++){
523 if(sinputFile[iF].Length()==0)continue;
524 // Draw the jet spectra
525 if(!inputList[iF])continue;
526 TH1F *hRec[nJets+1];
527 HistsFromSingleJets(cMaskRec,inputList[iF],nJets,hRec,iRebin);
528 /// TH1F *hRec[nJets] = (TH1F*)inputList[iF]->FindObject("fh1PtJetsRecIn"); hRec[nJets]->Rebin(iRebin);
529
530 hRec[nJets]->SetMarkerStyle(kFullCircle);
531 hRec[nJets]->SetXTitle("p_{T} (GeV/c)");
532 hRec[nJets]->SetMarkerColor(kColor[iF]);
533 hRec[nJets]->Scale(1./1.6); // deta
534 hRec[nJets]->Scale(1./hRec[nJets]->GetBinWidth(1));
535
536 hRec[nJets]->SetYTitle("dN/dp_{T}dy");
537
538 if(true){
539
540 ScaleNevent(hRec[nJets], inputFile[iF],AliAnalysisTaskJetServices::kSelected,0,iF); // <----- careful assumes iCen == iF
541 hRec[nJets]->SetYTitle("1/N_{evt} dN/dp_{T}");
542 // Printf("%d Ncall %f",iF,nColl[iF]);
543 // hRec[nJets]->Scale(1./nColl[iF]);
544 // hRec[nJets]->SetYTitle("1/(N_{coll}*(N_{evt}) dN/dp_{T}dy");
545 yMin = 1E-8;
546 yMax = 1;
547 }
548 if(sCen[iF].Length()){
549 leg1->AddEntry(hRec[nJets],Form("%s",sCen[iF].Data()),"P");
550 }
551
552 hJets[iF] = hRec[nJets];
553 c1->cd();
554
555
556 if(bLogLog)hRec[nJets]->SetAxisRange(5,100);
557 else hRec[nJets]->SetAxisRange(1,100);
558
559 if(iF==0){
560 if(yMax>0) hRec[nJets]->SetMaximum(yMax);
561 if(yMin>0) hRec[nJets]->SetMinimum(yMin);
562 hRec[nJets]->DrawCopy("P");
563 }
564 else{
565 if(sinputJet[iF].Length())hRec[nJets]->DrawCopy("Psame");
566 }
567 c1->Update();
568 if(getchar()=='q')return 1;
569 }
570
571
572 c1->Update();
573 TString picName = "jetspectrumPbPb";
574 if(bLogLog)picName += "LogLog";
575 leg1->Draw("same");
576 TLatex *txt = 0;
577 txt = new TLatex(5,0.1,"LHC2010 Pb+Pb #sqrt{s_{NN}} = 2.76 TeV");
578 txt->SetTextSize(gStyle->GetTextSize()*0.8);
579 txt->Draw();
580 ALICEWorkInProgress(c1,"02/15/2011");
581 c1->SaveAs(Form(cPrintMask,picName.Data()));
582 if(getchar()=='q')return 1;
583
584
585 c2->cd();
586 leg1->Draw("same");
587
588
589
590
591 picName = "jetspectrumLabels";
592 c2->SaveAs(Form(cPrintMask,picName.Data()));
593 if(getchar()=='q')return 1;
594
595 // scale with nColl
596 for(int iF = 0; iF < kMaxFiles;iF++){
597 if(sinputFile[iF].Length()==0)continue;
598 if(!hJets[iF])continue;
599 Printf("%d Ncoll %f",iF,nColl[iF]);
600 hJets[iF]->Scale(1./nColl[iF]);
601 hJets[iF]->SetYTitle("1/(N_{coll}*(N_{evt}) dN/dp_{T}dy");
602 }
603
604 const Int_t nPSmear = 9;
605 TH1F *hPythia[nPSmear];
606 // Fetch the smeared pythia spectra
607
608 TString pName;
609 c2->SetLogy();
610 TH1F* hRatio[kMaxFiles];
611
612 for(int i = 0;i < nPSmear;i++){
613 if(i==nPSmear-1){
614 pName = "h1JetPtCh";
615 hPythia[i] = GetPythia8Spectrum(pName.Data(),"~/alice/sim/pythia8/examples/output/LHCJets_2.75TeV_allpt_R04.root",1.0,iRebin);
616 }
617 else {
618 pName = Form("h1JetPtCh_C%02d_j%%d",i);
619 hPythia[i] = GetPythia8Spectrum(pName.Data(),"~/alice/sim/pythia8/examples/output/LHCJets_2.75TeV_allpt_R04.root",1.0,iRebin,nJets);
620 }
621 hPythia[i]->SetMarkerStyle(kOpenCircle);
622 hPythia[i]->SetMarkerColor(kBlack);
623 hPythia[i]->SetYTitle("1/(N_{coll}*N_{evt}) dN/dp_{T}dy");
624 hPythia[i]->SetAxisRange(0,100);
625 c1->cd();
626 hPythia[i]->Draw("CP");
627 hPythia[i]->SetMaximum(1E-02);
628 hPythia[i]->SetMinimum(1E-09);
629 for(int iF = 0; iF < kMaxFiles;iF++){
630 if(sinputFile[iF].Length()==0)continue;
631 if(!hJets[iF])continue;
632
633 c1->cd();
634 hJets[iF]->Draw("psame");
635
636 c2->cd();
637 hRatio[iF] = (TH1F*)hJets[iF]->Clone(Form("hRatio%d",iF));
638 hRatio[iF]->SetMaximum(100);
639 hRatio[iF]->SetMinimum(0.01);
640 hRatio[iF]->Divide(hPythia[i]);
641 hRatio[iF]->SetYTitle("Ratio");
642 if(iF==0) hRatio[iF]->DrawCopy("p");
643 else hRatio[iF]->DrawCopy("psame");
644 c1->Update();
645 c2->Update();
646 }
647
648 if(getchar()=='q')return 1;
649 picName = Form("jetspectrumRatioPbPb_Smear%d",i);
650 c2->SaveAs(Form(cPrintMask,picName.Data()));
651 picName = Form("jetspectrumPbPb_Smear%d",i);
652 c1->SaveAs(Form(cPrintMask,picName.Data()));
653 if(getchar()=='q')return 1;
654
655 }
656 c1->cd();
657 hPythia[nPSmear-1]->DrawCopy("P");
658 for(int i = 0;i < nPSmear-1;i++){
659 c1->cd();
660 hPythia[i]->SetMarkerStyle(kFullCircle);
661 hPythia[i]->SetMarkerColor(kRed-9+i);
662 hPythia[i]->DrawCopy("psame");
663
664 TH1F *hRatioS = (TH1F*)hPythia[i]->Clone(Form("hPythia_%d",i));
665 hRatioS->Divide(hPythia[nPSmear-1]);
666 hRatioS->SetMaximum(1E5);
667 hRatioS->SetMinimum(0.1);
668 hRatioS->SetYTitle("Ratio");
669 c2->cd();
670 if(i==0)hRatioS->DrawCopy("p");
671 else hRatioS->DrawCopy("psame");
672 c1->Update();
673 c2->Update();
674 if(getchar()=='q')return 1;
675 }
676
677 if(getchar()=='q')return 1;
678 picName = Form("jetSmearRatio");
679 c2->SaveAs(Form(cPrintMask,picName.Data()));
680 picName = Form("jetSmear");
681 c1->SaveAs(Form(cPrintMask,picName.Data()));
682 if(getchar()=='q')return 1;
683
684 // fetch the jet spectrum first
685
686
687 c1->cd();
688 TH1F *hPythiaRef = GetPythia8Spectrum(pName.Data(),"~/alice/sim/pythia8/examples/output/LHCJets_2.75TeV_allpt_R04.root",1.0,1);
689 hPythiaRef->Draw();
690 c1->Update();
691 Float_t fMinPt[kMaxFiles];
692 for(int iF = 0; iF < kMaxFiles;iF++){
693 if(nColl[iF]<=0)continue;
694 for(int i = hPythiaRef->GetNbinsX();i >0;i--){
695 if(hPythiaRef->GetBinContent(i)>1./nColl[iF]){
696 fMinPt[iF] = hPythiaRef->GetBinCenter(i);
697 break;
698 }
699 }
700 Printf("Ncoll %4.1f Min Pt %3.1f",nColl[iF],fMinPt[iF]);
701 }
702
703 char *cFile = "~/alice/sim/pythia8/examples/output/LHCJets_2.75TeV_allpt_R04.root";
704 TFile *fIn = (TFile*)gROOT->GetListOfFiles()->FindObject(cFile);
705 if(!fIn) fIn = TFile::Open(cFile);
706 if(!fIn)return 0;
707
708 TH1D* hPythia2[nPSmear];
709 for(int i = 0;i < nPSmear-1;i++){
710
711 TH2F *hCorrelation = 0;
712 for(int ij = 0;ij<2;ij++){
713 TH2F* hTmp = (TH2F*)fIn->Get(Form("h2SmearCorrelationCh_C%02d_j%d",i,ij));
714 if(hTmp&&!hCorrelation)hCorrelation = (TH2F*)hTmp->Clone(Form("%s_%d",hTmp->GetName(),i));
715 else if (hTmp) hCorrelation->Add(hTmp);
716 }
717 c1->cd();
718
719 hPythia2[i] = hCorrelation->ProjectionY(Form("hPythia2_%d",i),10,300);
720 hPythia2[i]->SetMarkerStyle(kFullCircle);
721 hPythia2[i]->SetMarkerColor(kRed-9+i);
722 hPythia2[i]->Rebin(iRebin);
723 hPythia2[i]->SetAxisRange(0,100);
724 hPythia2[i]->Scale(1./hPythia2[i]->GetBinWidth(1));
725 // hist->Scale(1./xsec);
726 hPythia2[i]->Scale(1./1.); // deta
727
728 // Double_t xtotal = 71.39; // xtotal for 7 TeV!!
729 Double_t xtotal = 62.01; // xtotal for 2.76 TeV!!
730 // scale with fraction of total NSD cross section
731 hPythia2[i]->Scale(1/xtotal);
732
733
734
735 c1->cd();
736 if(i==0)hPythia2[i]->DrawCopy("p");
737 else hPythia2[i]->DrawCopy("psame");
738
739 TH1F *hRatioS = (TH1F*)hPythia2[i]->Clone(Form("hPythia2_%d",i));
740 hRatioS->Divide(hPythia[nPSmear-1]);
741 hRatioS->SetMaximum(1E5);
742 hRatioS->SetMinimum(0.1);
743 hRatioS->SetYTitle("Ratio");
744 c2->cd();
745 if(i==0)hRatioS->DrawCopy("p");
746 else hRatioS->DrawCopy("psame");
747 c1->Update();
748 c2->Update();
749 if(getchar()=='q')return 1;
750 }
751 return 0;
752}
753
754
755
756void set_plot_style() {
757 const Int_t NRGBs = 5;
758 const Int_t NCont = 255;
759
760 Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
761 Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
762 Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
763 Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
764 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
765 gStyle->SetNumberContours(NCont);
766}
767
768
769Int_t PlotJetBFluctuations(){
770 // plot the diffent background estimates
771
772 const int nCen = 4;
773 const Float_t fCentLo[nCen] = {0,10,30,50};
774 const Float_t fCentUp[nCen] = {10,30,50,80};
775 TH2F *hFrame = new TH2F("hFrame",";#delta p_{T} (GeV/c);Probability/GeV",200,-70,70,100,1E-5,5);
776 hFrame->SetTitleOffset(1.5,"Y");
777 hFrame->SetTitleOffset(1.5,"X");
778 hFrame->SetLabelSize(hFrame->GetLabelSize("Y")*0.9,"Y");
779 hFrame->SetLabelSize(hFrame->GetLabelSize("X")*0.9,"X");
780 hFrame->SetTitleSize(hFrame->GetTitleSize("Y")*0.7,"Y");
781 hFrame->SetTitleSize(hFrame->GetTitleSize("X")*0.7,"X");
782
783
784 TString printType = "png";
785 TString tmpName;
786
787 TCanvas *c1 = new TCanvas("c11","c1",600,600);
788 c1->SetLogy();
789
790 TFile::SetCacheFileDir("/tmp/");
791
792 // Martha, single particle jets
793 TFile *fM = TFile::Open("http://qgp.uni-muenster.de/~stevero/tmp/jets/PWG4_JetTasksOutput_AOD_EmbeddingSingleTrack.root","CACHEREAD");
794 TH1D *hDeltaPtM[nCen] = {0};
795
796 // select
797 Float_t fMinPtM = 40;
798 Float_t fMaxPtM = 80;
799 int iB = 2;
800
801 /*
802 0: 0-10%
803 1: 10-30%
804 2: 30-50%
805 3: 50-80%
806 */
807
808
809
810
811 for(int ic = 0;ic < nCen;ic++){
812 tmpName = Form("PWG4_BkgFluctCent%dB%d",ic,iB);
813 TDirectory *dir = (TDirectory*)fM->Get(tmpName.Data());
814 if(!dir)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
815 tmpName = Form("taskBkgFluctCent%dB%d",ic,iB);
816 TList *list = (TList*)dir->Get(tmpName.Data());
817 if(!list)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
818
819 TH3F *h3Tmp = (TH3F*)list->FindObject("fh3PtDeltaPtArea");
820 hDeltaPtM[ic] = h3Tmp->ProjectionY(Form("hDeltaM%d",ic),h3Tmp->GetXaxis()->FindBin(fMinPtM),h3Tmp->GetXaxis()->FindBin(fMaxPtM),
821 0,-1,"E");
822 hDeltaPtM[ic]->Rebin(10);
823 Float_t fScale = hDeltaPtM[ic]->Integral("width");
824 if(fScale)hDeltaPtM[ic]->Scale(1./fScale);
825 hDeltaPtM[ic]->SetMarkerStyle(kFullSquare);
826 hDeltaPtM[ic]->SetMarkerColor(kGreen+2);
827 hDeltaPtM[ic]->SetLineColor( hDeltaPtM[ic]->GetMarkerColor());
828 }
829
830 // fetch the BiAs
831
832 TH1D *hBiaL[nCen];
833
834 TFile *fL = TFile::Open("~/alice/jets/macros/corrections/tmp/pwg4plots.root");
835
836 for(int ic = 0;ic < nCen;ic++){
837 if(iB==1)tmpName = "BiA sa antikt centrality";
838 else if (iB==2)tmpName = "BiA va antikt centrality";
839 TH2F *h2Tmp = (TH2F*)fL->Get(tmpName.Data());
840 if(!h2Tmp)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
841 Int_t ibLo = h2Tmp->GetXaxis()->FindBin(fCentLo[ic]);
842 Int_t ibUp = h2Tmp->GetXaxis()->FindBin(fCentUp[ic])-1;
843 Printf("Line:%d bin %d - %d",__LINE__,ibLo,ibUp);
844 hBiaL[ic] = h2Tmp->ProjectionY(Form("hBiaL%d",ic),ibLo,ibUp,"E");
845 Float_t fScale = hBiaL[ic]->Integral("width");
846 if(fScale) hBiaL[ic]->Scale(1./fScale);
847 hBiaL[ic]->SetMarkerStyle(kFullCircle);
848 hBiaL[ic]->SetMarkerColor(kBlue);
849 hBiaL[ic]->SetLineColor( hBiaL[ic]->GetMarkerColor());
850 }
851
852 TH1D *hBiaRC[nCen];
853
854 for(int ic = 0;ic < nCen;ic++){
855 if(iB==1)tmpName = "BiA sa RC centrality";
856 else if (iB==2)tmpName = "BiA va RC centrality";
857 TH2F *h2Tmp = (TH2F*)fL->Get(tmpName.Data());
858 if(!h2Tmp)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
859 Int_t ibLo = h2Tmp->GetXaxis()->FindBin(fCentLo[ic]);
860 Int_t ibUp = h2Tmp->GetXaxis()->FindBin(fCentUp[ic])-1;
861 Printf("Line:%d bin %d - %d",__LINE__,ibLo,ibUp);
862 hBiaRC[ic] = h2Tmp->ProjectionY(Form("hBiaRC%d",ic),ibLo,ibUp,"E");
863 Float_t fScale = hBiaRC[ic]->Integral("width");
864 if(fScale) hBiaRC[ic]->Scale(1./fScale);
865 hBiaRC[ic]->SetMarkerStyle(kFullCircle);
866 hBiaRC[ic]->SetMarkerColor(kRed);
867 hBiaRC[ic]->SetLineColor( hBiaRC[ic]->GetMarkerColor());
868 }
869
870
871
872 c1->SetLogy(0);
873 c1->SetLogz();
874 TLatex *txt = new TLatex();
875 txt->SetTextFont(gStyle->GetTextFont());
876 txt->SetTextSize(gStyle->GetTextSize()*0.6);
877
878 TLatex *txt2 = new TLatex();
879 txt2->SetTextFont(gStyle->GetTextFont());
880 txt2->SetTextSize(gStyle->GetTextSize()*0.7);
881 txt2->SetTextAlign(22);
882 txt2->SetTextColor(kRed);
883
884
885 if(iB==1)tmpName = "background sa vs multiplicity";
886 else if(iB==2)tmpName = "background va vs multiplicity";
887 TH2F *h2RhoVsMult = (TH2F*)fL->Get(tmpName.Data());
888 if(!h2RhoVsMult)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
889 c1->SetMargin(0.15,0.18,0.2,0.10);
890
891 h2RhoVsMult->SetTitleOffset(1.5,"Y");
892 h2RhoVsMult->SetTitleOffset(1.5,"X");
893 h2RhoVsMult->SetLabelSize(h2RhoVsMult->GetLabelSize("Y")*0.9,"Y");
894 h2RhoVsMult->SetLabelSize(h2RhoVsMult->GetLabelSize("X")*0.9,"X");
895 h2RhoVsMult->SetTitleSize(h2RhoVsMult->GetTitleSize("Y")*0.7,"Y");
896 h2RhoVsMult->SetTitleSize(h2RhoVsMult->GetTitleSize("X")*0.7,"X");
897 h2RhoVsMult->SetXTitle("input tracks");
898 h2RhoVsMult->SetYTitle("#rho (GeV/unit area)");
899 h2RhoVsMult->SetAxisRange(0,3000.);
900 h2RhoVsMult->Draw("colz");
901 txt->Draw();
902 txt->DrawLatex(100,180,"LHC 2010 Pb+Pb Run #sqrt{s_{NN}} = 2.76 TeV");
903 txt2->DrawLatex(800,150,"ALICE Performance");
904 txt2->DrawLatex(800,140,"01/03/2011");
905 c1->Update();
906 c1->SaveAs(Form("rhovsmult_B%d.%s",iB,printType.Data()));
907 if(getchar()=='q')return 1;
908
909 if(iB==1)tmpName = "background sa vs centrality";
910 else if(iB==2)tmpName = "background va vs centrality";
911 TH2F *h2RhoVsCent = (TH2F*)fL->Get(tmpName.Data());
912 if(!h2RhoVsCent)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
913 h2RhoVsCent->SetXTitle("centrality (%)");
914 h2RhoVsCent->SetYTitle("#rho (GeV/unit area)");
915
916 h2RhoVsCent->SetTitleOffset(1.5,"Y");
917 h2RhoVsCent->SetTitleOffset(1.5,"X");
918 h2RhoVsCent->SetLabelSize(h2RhoVsCent->GetLabelSize("Y")*0.9,"Y");
919 h2RhoVsCent->SetLabelSize(h2RhoVsCent->GetLabelSize("X")*0.9,"X");
920 h2RhoVsCent->SetTitleSize(h2RhoVsCent->GetTitleSize("Y")*0.7,"Y");
921 h2RhoVsCent->SetTitleSize(h2RhoVsCent->GetTitleSize("X")*0.7,"X");
922 h2RhoVsCent->SetAxisRange(3,200,"Y");
923 h2RhoVsCent->Draw("colz");
924 txt->DrawLatex(20,180,"LHC 2010 Pb+Pb Run #sqrt{s_{NN}} = 2.76 TeV");
925 txt2->DrawLatex(50,150,"ALICE Performance");
926 txt2->DrawLatex(50,140,"01/03/2011");
927 c1->Update();
928 c1->SaveAs(Form("rhovscent_B%d.%s",iB,printType.Data()));
929 if(getchar()=='q')return 1;
930
931 // fetch the data from bastian...
932 Float_t fMinPtB = fMinPtM;
933 Float_t fMaxPtB = fMaxPtM;
934
935
936 // the embbedded jets, carefull, take only above 60 GeV
937 // c = all jets, d = leading jets, e = single tracks
938 TFile *fB1 = TFile::Open("~/alice/jets/macros/corrections/tmp/PWG4_JetTasksOutput_e.root");
939 TH1D *hDeltaPtB1[nCen] = {0};
940 for(int ic = 0;ic < nCen;ic++){
941 tmpName = Form("PWG4_JetResponse_%s_%s%02d_B%d_Filter00256_Cut%05d_Skip00","clusters", "ANTIKT", (Int_t)(0.4*10), iB, (Int_t)(0.15*1000));
942 TDirectoryFile *df = dynamic_cast<TDirectoryFile*> (fB1->Get(tmpName.Data()));
943 if(!df)Printf("%d %s not found",__LINE__,tmpName.Data());
944 tmpName.ReplaceAll("PWG4_JetResponse","jetresponse");
945 TList *list = dynamic_cast<TList*> (df->Get(tmpName.Data()));
946 if(!list)Printf("%d %s not found",__LINE__,tmpName.Data());
947 tmpName = Form("pt_smearing%d",ic+1);
948 TH2F *hTmp = (TH2F*)list->FindObject(tmpName.Data());
949 if(!hTmp)Printf("%d %s not found",__LINE__,tmpName.Data());
950 int ibLo = hTmp->GetYaxis()->FindBin(fMinPtB);
951 int ibUp = hTmp->GetYaxis()->FindBin(fMaxPtB)-1;
952 hDeltaPtB1[ic] = hTmp->ProjectionX(Form("fHistDeltaPtB1_c%d",ic),ibLo,ibUp,"E");
953 hDeltaPtB1[ic]->SetMarkerStyle(kFullSquare);
954 hDeltaPtB1[ic]->SetMarkerColor(kBlue);
955 hDeltaPtB1[ic]->SetLineColor(hDeltaPtB1[ic]->GetMarkerColor());
956 hDeltaPtB1[ic]->Rebin(2);
957 Float_t fScale = hDeltaPtB1[ic]->Integral("width");
958 if(fScale) hDeltaPtB1[ic]->Scale(1./fScale);
959 }
960
961
962 c1->SetLogy();
963 c1->SetMargin(0.15,0.05,0.2,0.05);
964
965 TF1 *gaus = new TF1("gaus","gaus",-60,2);
966 TF1 *gaus2 = new TF1("gaus2","gaus",-60,2);
967 for(int ic = 0;ic < nCen;ic++){
968 TLegend *leg1 = new TLegend(0.2,0.78,0.3,0.93);
969 leg1->SetHeader(Form("Pb+Pb %2.0f-%2.0f%% R = 0.4 (B%d)",fCentLo[ic],fCentUp[ic],iB));
970 leg1->SetFillColor(0);
971 leg1->SetTextFont(gStyle->GetTextFont());
972 leg1->SetTextSize(gStyle->GetTextSize()*0.6);
973 leg1->SetBorderSize(0);
974 hFrame->DrawCopy();
975
976 /*
977 hBiaL[ic]->DrawCopy("psame");
978 leg1->AddEntry(hBiaL[ic],Form("BiA anti-k_{T}"),"P");
979 */
980
981 hBiaRC[ic]->DrawCopy("psame");
982 gaus->SetLineColor(hBiaRC[ic]->GetLineColor());
983 hBiaRC[ic]->Fit(gaus,"R0");
984 gaus->SetRange(-40.,40.);
985 gaus->Draw("same");
986 leg1->AddEntry(hBiaRC[ic],Form("BiA random cones (excl. leading jets)"),"P");
987 leg1->AddEntry(gaus,Form("LHS Gaus fit: #mu = %2.2f, #sigma = %2.2f",gaus->GetParameter(1),gaus->GetParameter(2)),"L");
988
989
990
991 // hDeltaPtB1[ic]->Draw("psame");
992
993
994 hDeltaPtM[ic]->DrawCopy("psame");
995 gaus2->SetLineColor(hDeltaPtM[ic]->GetLineColor());
996 hDeltaPtM[ic]->Fit(gaus2,"R0");
997 gaus2->SetRange(-40.,40.);
998 // gaus2->Draw("same");
999 leg1->AddEntry(hDeltaPtM[ic],Form("anti-k_{T} single tracks %2.0f-%2.0f GeV",fMinPtM,fMaxPtM),"P");
1000 leg1->Draw();
1001 txt2->SetNDC();
1002 txt2->DrawLatex(0.5,0.97,"ALICE Performance 01/03/2011");
1003 c1->Update();
1004 c1->SaveAs(Form("deltaPt_B%d_cen%02d.%s",iB,ic,printType.Data()));
1005 if(getchar()=='q')return 1;
1006
1007
1008 }
1009
1010
1011
1012 return 0;
1013
1014
1015}
1016
1017void ScaleH1(TH1* h,char* cX,char* cY,Float_t fScale,Bool_t bWidth){
1018 if(!h)return;
1019 if(!fScale)return;
1020
1021 h->Scale(fScale);
1022 if(bWidth){
1023 for(int ib = 1;ib <= h->GetNbinsX();ib++){
1024 Float_t val = h->GetBinContent(ib);
1025 Float_t err = h->GetBinError(ib);
1026 Float_t w = h->GetBinWidth(ib);
1027 h->SetBinContent(ib,val/w);
1028 h->SetBinError(ib,err/w);
1029 // Printf("width %f",w);
1030 }
1031 }
1032 h->SetXTitle(cX);
1033 h->SetYTitle(cY);
1034
1035}
1036
1037
1038void ScaleH2(TH2* h,char* cX,char* cY,char* cZ,Float_t fScale,Bool_t bWidth){
1039 if(!h)return;
1040 if(!fScale)return;
1041
1042 h->Scale(fScale);
1043 if(bWidth){
1044 for(int ibx = 1;ibx <= h->GetNbinsX();ibx++){
1045 for(int iby = 1;iby <= h->GetNbinsY();iby++){
1046 Float_t val = h->GetBinContent(ibx,iby);
1047 Float_t err = h->GetBinError(ibx,iby);
1048 Float_t wx = h->GetXaxis()->GetBinWidth(ibx);
1049 Float_t wy = h->GetYaxis()->GetBinWidth(iby);
1050 h->SetBinContent(ibx,iby,val/(wx*wy));
1051 h->SetBinError(ibx,iby,err/(wx*wy));
1052 // Printf("width %f",w);
1053 }
1054 }
1055 }
1056 h->SetXTitle(cX);
1057 h->SetYTitle(cY);
1058 h->SetZTitle(cZ);
1059
1060}
1061
1062void SetStyleH1(TH1 *h,Int_t color,Int_t fillStyle,Int_t markerStyle){
1063 if(color){
1064 h->SetLineColor(color);
1065 h->SetLineWidth(3);
1066 h->SetMarkerColor(color);
1067 h->SetFillColor(color);
1068 }
1069 if(fillStyle)h->SetFillStyle(fillStyle);
1070 if(markerStyle)h->SetMarkerStyle(markerStyle);
1071}
1072
1073void HistsFromSingleJets(char* cMask,TList *list,Int_t nJets,TH1F** h,Int_t iRebin,int iFirst){
1074
1075 for(int ij = iFirst;ij < nJets;ij++){
1076 h[ij] = (TH1F*)list->FindObject(Form(cMask,ij));
1077
1078 if(h[ij]){
1079 if(iRebin>1){
1080 h[ij] = (TH1F*)h[ij]->Clone();
1081 h[ij]->SetDirectory(gROOT);
1082 h[ij]->Rebin(iRebin);
1083 }
1084 if(ij==iFirst){
1085 h[nJets] = (TH1F*)h[ij]->Clone(Form(cMask,nJets));
1086 }
1087 else{
1088 h[nJets]->Add(h[ij]);
1089 }
1090 }
1091 else{
1092 Printf("%s not found",Form(cMask,ij));
1093 }
1094 }
1095
1096}
1097
1098
1099void HistsFromSingleJets(char* cMask,TList *list,Int_t nJets,TH2F** h,Int_t iRebin,int iFirst){
1100
1101 for(int ij = iFirst;ij < nJets;ij++){
1102 h[ij] = (TH2F*)list->FindObject(Form(cMask,ij));
1103 if(h[ij]){
1104 if(iRebin>1){
1105 h[ij] = (TH2F*)h[ij]->Clone();
1106 h[ij]->SetDirectory(gROOT);
1107 h[ij]->Rebin(iRebin);
1108 }
1109 if(ij==iFirst){
1110 h[nJets] = (TH2F*)h[ij]->Clone(Form(cMask,nJets));
1111 }
1112 else{
1113 h[nJets]->Add(h[ij]);
1114 }
1115 }
1116 else{
1117 Printf("%s not found",Form(cMask,ij));
1118 }
1119 }
1120
1121}
1122
1123void ScaleNevent(TH1* h,TFile *fIn,Int_t iCond,Int_t iMC,Int_t iCen){
1124 // fetch the trigger histos
1125
1126 TDirectory *dir = (TDirectory*)fIn->Get("PWG4_services");
1127 TList *list1 = (TList*)dir->Get("pwg4serv");
1128
1129 TH2F* hTriggerCount =(TH2F*)list1->FindObject("fh2TriggerCount");
1130 Float_t nevents = hTriggerCount->GetBinContent(hTriggerCount->FindBin(0+iCen,iCond)); // we have to scale with number of triggers?
1131
1132 Float_t xsection = 1;
1133 if(iMC==7000){
1134 // 7 TeV sigma ND = 48.85 sigma NEL = 71.39 (Pythia 8... confirm with 6
1135 xsection = 71.39;
1136 Printf("MC Scaling setting nevents=%f to 1, xsection = %E",nevents,xsection);
1137 nevents = 1;
1138 }
1139
1140 Printf("Scaling %s with number of event %f",h->GetName(),nevents);
1141 if(nevents){
1142 Float_t scalef = 1./nevents/xsection;
1143 Printf("scale factor %E",scalef);
1144 h->Scale(scalef);
1145 }
1146
1147}
1148
814a52a0 1149TF1* FitLHSgaus(TH1D *hDeltaPt, double &sigma, double &sigma_error, double &mean)
1150{
1151 hDeltaPt->Fit("gaus",-50.,0.);
1152 TF1 *f1 = hDeltaPt->GetFunction("gaus");
1153 sigma = f1->GetParameter(2);
1154 mean = f1->GetParameter(1);
1155 hDeltaPt->Fit("gaus","""","""",mean-3.*sigma,0.);
1156 f1 = hDeltaPt->GetFunction("gaus");
1157 mean = f1->GetParameter(1);
1158 sigma = f1->GetParameter(2);
1159 sigma_error = f1->GetParError(2);
1160
1161 return f1;
1162}