]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/macros/jetspectrum/PlotNote.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fix bug from last commit - Marta V.
[u/mrichter/AliRoot.git] / PWGJE / macros / jetspectrum / PlotNote.C
Content-type: text/html ]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/macros/jetspectrum/PlotNote.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:


500 - Internal Server Error
Malformed UTF-8 character (fatal) at (eval 5) line 1, <$fd> line 6676.
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"
1354c5f3 34#include "utils.C"
e2081c22 35
36using namespace std;
37
38Int_t PlotSpectraPbPb();
5f5e3e77 39Int_t PlotSpectraPbPb2(UInt_t iPlotFlag = 0xFF,UInt_t iCenFlag = 0xFF,UInt_t iSpecFlag = 0xFF,Int_t iALICEType =100);
e2081c22 40Int_t PlotJetBFluctuations();
1354c5f3 41Int_t PlotJetBFluctuations2(UInt_t iPlotFlag = 0xFFFFF,UInt_t iPlotType = 0xFFFFF,Float_t inPtLo = 50,Float_t inPtUp = 100,Int_t iALICEType = 100);
5f5e3e77 42Int_t PlotSubtract(Int_t iALICEtype = 100);
e2081c22 43Int_t PlotJetQA();
5f5e3e77 44TH1* CorrectForEff(TH1 *hTrack);
e2081c22 45// helpers
46
8dc333da 47Int_t fDebug = 0;
e2081c22 48
49enum {kMaxJets = 4};
50void HistsFromSingleJets(char* cMask,TList *list,Int_t nJets,TH1F** h,Int_t iRebin = 1,int iFirst = 0);
51void HistsFromSingleJets(char* cMask,TList *list,Int_t nJets,TH2F** h,Int_t iRebin = 1,int iFirst = 0);
52void ScaleNevent(TH1 *h,TFile *fIn,Int_t iCond = 5,Int_t iMC = 0,Int_t iCen = 0);
53void ScaleH1(TH1* h,char* cX = "",char* cY = "",Float_t fScale = 1,Bool_t bWidth = true);
54void SetStyleH1(TH1 *h,Int_t color = 0,Int_t fillStyle = 0,Int_t markerStyle = 0);
55void ScaleH2(TH2* h,char* cX = "",char* cY = "",char* cZ = "",Float_t fScale = 1,Bool_t bWidth = true);
56TH1F* GetPythia8Spectrum(const char *cHist,const char *cPythia,Float_t deta,Int_t iRebin = 1,Int_t iMask = 0);
8dc333da 57Double_t GetPoissonFluctuation(TH1 *h1,Double_t areaIn,Double_t areaJet);
8f206bf0 58TObject* GetObjectFromFile(const char *cName,const char *cFile,const char* cDir,const char* cRep2 = "pwg4",const char* cRep1 = "PWG4_");
59TF1* FitLHSgaus(TH1D *hDeltaPt, double minPt = -60., double maxPt = 5., int minIterations = 3, int maxIterations = 20, double maxDeltaMean = 0.01, int mode=1, double minBoundSigma = 3., double maxBoundSigma = 0.5 );
c08c3ad2 60Double_t Gamma(Double_t *x,Double_t *par);
61Double_t Gamma0(Double_t *x,Double_t *par);
8dc333da 62TH2F* GetTH2PlotB(const char *cPath,Int_t embType=0, Int_t classType=0, Int_t cl=-1, Int_t rp=-1);
8f206bf0 63void SetHistoAttributes(TH1* h1,Int_t iMarker = kFullCircle,Int_t iColor = kBlack,Float_t fMarkerSize = 1.7);
64void SetHistAxisStyle(TH1* h1,Float_t titleOffsetX = 0,Float_t titleSizeX= 0,Float_t labelSizeX= 0,
65 Float_t titleOffsetY = 0,Float_t titleSizeY = 0,Float_t labelSizeY = 0);
8dc333da 66void SetGraphAttributes(TGraph* gr,Int_t iMarker = kFullCircle,Int_t iColor = kBlack);
c08c3ad2 67
8dc333da 68void CloseFiles();
8f206bf0 69void DrawDate();
5f5e3e77 70void DrawALICE(TCanvas *c,Int_t iType = 11,Float_t xCenter = 0.90,Float_t yCenter = 0.9,TString cExtra2 = "",bool bExtraBox = kTRUE);
e2081c22 71const char *cPrintMask = "110116_%s.eps";
72void set_plot_style();
73TCanvas *cTmp = 0;
74
8f206bf0 75TString picPrefix = "notePlot_";
76TString picSuffix = "eps";
8dc333da 77TString picName = "";
78
8f206bf0 79
80
81
e2081c22 82TH1F* GetPythia8Spectrum(const char *cHist,const char *cPythia,Float_t deta,Int_t iRebin,Int_t iMask){
83 TDirectory *opwd = gDirectory;
84 TFile *fPythia = (TFile*)gROOT->GetListOfFiles()->FindObject(cPythia);
85
86 if(!fPythia) fPythia = TFile::Open(cPythia);
87 opwd->cd();
88
89 static int iCount = 0;
90 TH1F* hist = 0;
91 if(iMask<=0){
92 hist = (TH1F*)fPythia->Get(cHist);
93 if(hist)hist = (TH1F*)hist->Clone(Form("%s_%d",hist->GetName(),iCount));
94 }
95 else{
96 for(int i = 0;i < iMask;i++){
97 TH1F *hTmp = (TH1F*)fPythia->Get(Form(cHist,i));
98 if(!hTmp){
99 Printf("%s not found",Form(cHist,i));
100 return 0;
101 }
102 if(!hist)hist = (TH1F*)hTmp->Clone(Form("%s_%d_%d",hTmp->GetName(),iMask,iCount));
103 else hist->Add(hTmp);
104 }
105 }
106
107 if(!hist){
108 Printf("%s not found",cHist);
109 return 0;
110 }
111 // fetch the cross section
112 TH1F *hxsec = (TH1F*)fPythia->Get("h1Xsec");
113
114 Float_t xsec = hxsec->GetBinContent(1);
115 Printf("%d xsec = %1.3E",__LINE__,xsec);
116 // if(xsec==0)xsec = 40.79; // tmp hack
117 hist->Rebin(iRebin);
118 hist->Scale(1./hist->GetBinWidth(1));
119 // hist->Scale(1./xsec);
120 hist->Scale(1./deta);
121
122 // Double_t xtotal = 71.39; // xtotal for 7 TeV!!
123 Double_t xtotal = 62.01; // xtotal for 7 TeV!!
124
125 // scale with fraction of total NSD cross section
126 hist->Scale(1/xtotal);
127 Printf("%d fraction = %1.3E of total cross section %1.3E",__LINE__,xsec/xtotal,xtotal);
128 hist->SetXTitle("p_{T} (GeV/c)");
129 hist->SetYTitle("1/N dN/dp_{T}dy");
130 hist->SetDirectory(opwd);
131 iCount++;
132 return hist;
133
134}
135
c08c3ad2 136
8dc333da 137
138
e2081c22 139void PlotNote(){
8dc333da 140 set_plot_style();
8dc333da 141 // QM Plots
142
8f206bf0 143 // BD
5f5e3e77 144 PlotSubtract(12); // do not call before others changes style
8dc333da 145
5f5e3e77 146 // BF 01a-d Random cones
1354c5f3 147 // PlotJetBFluctuations2((1<<0)|(1<<1)|(1<<2),1<<0,50,100,13);
8dc333da 148
5f5e3e77 149 // BF 02 a-d jet, take out unquenched
150 // PlotJetBFluctuations2((1<<4)|(1<<6),1<<0);
8dc333da 151
5f5e3e77 152 // BF 03 a-d
1354c5f3 153 // PlotJetBFluctuations2((1<<3),1<<0);
8dc333da 154
5f5e3e77 155 // BF 04 a-d 05 06
1354c5f3 156 // PlotJetBFluctuations2((1<<0)|(1<<3)|(1<<4)|(1<<7)|(1<<8),1<<0,50,100,13);
8dc333da 157
5f5e3e77 158 // BF 07 08 (vs. Mult)
1354c5f3 159 // PlotJetBFluctuations2((1<<0)|(1<<1)|(1<<3)|(1<<4)|(1<<7)|(1<<8),1<<1,50,100,13);
5f5e3e77 160
161 // BF08
162 // PlotJetBFluctuations2((1<<8)|(1<<11),1<<1);
8f206bf0 163
5f5e3e77 164 // BF 11a-d
1354c5f3 165 // PlotJetBFluctuations2((1<<3),1<<0|1<<2);
8f206bf0 166
5f5e3e77 167 // BF 12a-d Random cones vs RP
8f206bf0 168 // PlotJetBFluctuations2((1<<0),1<<0|1<<2);
8dc333da 169
5f5e3e77 170 // BF 13 + 14 may add random cones here
171 // PlotJetBFluctuations2((1<<3)|(1<<8),1<<1|1<<2);
172
173
174 // BF 13 + 15 alternative
175 // PlotJetBFluctuations2((1<<0)|(1<<1)|(1<<8),1<<1|1<<2);
176
177 // BF 19
178 // PlotJetBFluctuations2((1<<9),(1<<1)|(1<<2));
179
180
8dc333da 181
8f206bf0 182 // JS 01 + 02
5f5e3e77 183 // PlotSpectraPbPb2(1<<0|1<<1,1<<0|1<<1|1<<2|1<<3,0,12);
8dc333da 184
8f206bf0 185 // not for external usage with tracks...
5f5e3e77 186 // PlotSpectraPbPb2(1<<1,1<<0|1<<1|1<<2|1<<3,1);
8dc333da 187
e2081c22 188}
189
190TList *GetList(const char *cFile,const char *cName){
191 TDirectory *opwd = gDirectory;
192 TFile *fIn = (TFile*)gROOT->GetListOfFiles()->FindObject(cFile);
193 TList *list = 0;
194 if(!fIn) fIn = TFile::Open(cFile);
195 if(!fIn)return list;
196 opwd->cd();
197
198 TDirectory *dir = (TDirectory*)fIn->Get(Form("PWG4_%s",cName));
199 if(!dir){
200 Printf("GetList: Directory PWG4_%s not found",cName);
201 return list;
202 }
203 list = (TList*)dir->Get(Form("pwg4%s",cName));
204 if(!list){
205 Printf("GetList: list pwg4%s not found",cName);
206 return list;
207 }
208 return list;
209
210}
211
212Int_t PlotJetQA(){
213
214 set_plot_style();
215 gStyle->SetTitleYSize(0.7*gStyle->GetTitleYSize());
216 gStyle->SetTitleXSize(0.7*gStyle->GetTitleXSize());
217 gStyle->SetTitleXOffset(1.2*gStyle->GetTitleXOffset());
218 gStyle->SetTitleYOffset(1.5*gStyle->GetTitleYOffset());
219 gStyle->SetPadRightMargin(1.1*gStyle->GetPadRightMargin());
220 gStyle->SetPadBottomMargin(0.7*gStyle->GetPadBottomMargin());
221
222
223 Bool_t isPP = false;
224
225 TString printType = "eps";
226 const Int_t nCen = 5;
227 TString sCen[nCen] = {" 0-80"," 0-10%%","10-30","30-50","50-80"};
228 const Int_t nJets = 3;
229
230
231 TString cAdd = "Rec"; // decide whcih type to take Rec RecFull Gen GenFull
232 TString cAddFull = "RecFull"; // decide whcih type to take Rec RecFull Gen GenFull
233
234 TCanvas *c1 = new TCanvas();
235
236 TString sFile = "~/alice/data/analysis/train_maf/output_110216.root";
237 TString sinputName = "spec2_clustersAOD_ANTIKT04_B2_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000";
1354c5f3 238 TString sinputJet = "%% Pb-Pb Anti-k_{T} R = 0.4 (B2)";
e2081c22 239 TString sinputPrint = "PbPb_antikt04_B2";
1354c5f3 240 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";
e2081c22 241
242 if(isPP){
243 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";
244 sCen[0] = "";
245 sFile = "~/alice/jets/macros/batch/output.root";
246 }
247
248
249 TFile *fIn = TFile::Open(sFile.Data());
250
251
252 TDirectory *dir = (TDirectory*)fIn->Get("PWG4_services");
253 TList *list1 = (TList*)dir->Get("pwg4serv");
254
255 TH2F* hTriggerCount =(TH2F*)list1->FindObject("fh2TriggerCount");
256
257 for(int ic = 0;ic < (isPP?1:nCen);ic++){
258
259 TString cName = Form("PWG4_%s_Class%02d",sinputName.Data(),ic);
260 TDirectory *inputDir = (TDirectory*)fIn->Get(cName.Data());
261
262 Float_t fNevents = hTriggerCount->GetBinContent(hTriggerCount->FindBin(ic,AliAnalysisTaskJetServices::kSelected));
263
264 Printf(" %s: %10d events",sCen[ic].Data(),(Int_t)fNevents);
265
266 if(!inputDir){
267 Printf("Dir %s not found",cName.Data());
268 continue;
269 }
270 cName = Form("pwg4%s_Class%02d",sinputName.Data(),ic);
271 TList *inputList = (TList*)inputDir->Get(cName.Data());
272 if(!inputList){
273 Printf("List %s not found",cName.Data());
274 continue;
275 }
276
277 TH2F *h2EtaPtLead[nJets+2] = {0,};
278 TH2F *h2EtaPtLeadFull[nJets+2] = {0,};
279 TH2F *h2EtaPtAll[nJets+2] = {0,};
280 HistsFromSingleJets(Form("fh2EtaPt%s_j%%d",cAdd.Data()),inputList,1,h2EtaPtLead);
281 HistsFromSingleJets(Form("fh2EtaPt%s_j%%d",cAddFull.Data()),inputList,1,h2EtaPtLeadFull);
282 HistsFromSingleJets(Form("fh2EtaPt%s_j%%d",cAddFull.Data()),inputList,nJets+1,h2EtaPtAll,1,nJets);
283
284 c1->SetLogz();
285 for(int i = 0; i< nJets+2;i++)Printf("Lead > %d %p",i,h2EtaPtLead[i]);
286 for(int i = 0; i< nJets+2;i++)Printf("LeadFull > %d %p",i,h2EtaPtLeadFull[i]);
287 for(int i = 0; i< nJets+2;i++)Printf("All > %d %p",i,h2EtaPtAll[i]);
288 h2EtaPtLead[0]->Draw("colz");
289 c1->Modified();
290 c1->Update();
291 c1->SaveAs(Form("%s_%s_Cen%02d_EtaPtLead2D.%s",sinputPrint.Data(),cAdd.Data(),ic,printType.Data()));
292 if(!gROOT->IsBatch())
293 if(getchar()=='q')return 1;
294
295 Float_t fStep = 20;
296 Float_t fStepLo = 0;
297 Float_t fStepUp = fStepLo+fStep;
298 Float_t fMax = 160;
299
300 while(fStepUp<fMax){
301 Int_t iBinLo = h2EtaPtLead[0]->GetYaxis()->FindBin(fStepLo);
302 Int_t iBinUp = h2EtaPtLead[0]->GetYaxis()->FindBin(fStepUp)-1;
303 TH1D *hProjLead = h2EtaPtLead[0]->ProjectionX(Form("hEtaPtLead_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
304 TH1D *hProjLeadFull = h2EtaPtLeadFull[0]->ProjectionX(Form("hEtaPtLeadFull_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
305 TH1D *hProjAll = h2EtaPtAll[nJets+1]->ProjectionX(Form("hEtaPtAll_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
306
307 ScaleH1(hProjLead,"#eta","1/N_{evt} dN_{jet}/d#eta",1./fNevents,true);
308 SetStyleH1(hProjLead,kBlue,0,kFullCircle);
309 ScaleH1(hProjLeadFull,"#eta","1/N_{evt} dN_{jet}/d#eta",1./fNevents,true);
310 SetStyleH1(hProjLeadFull,kRed,0,kFullCircle);
311 ScaleH1(hProjAll,"#eta","1/N_{evt} dN_{jet}/d#eta",1./fNevents,true);
312 SetStyleH1(hProjAll,kGray,1001,kFullCircle);
313
314
315 hProjAll->DrawCopy("hist");
316 hProjLead->DrawCopy("histsame");
317 hProjLeadFull->DrawCopy("histsame");
318
319 Float_t ptUp = fStepUp-h2EtaPtLead[0]->GetYaxis()->GetBinWidth(iBinUp);
320
321 Printf("%3.0f-%3.0f",fStepLo,ptUp);
322 TLatex *txt = 0;
323 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));
324 txt->SetTextSize(gStyle->GetTextSize());
325 txt->Draw();
326 c1->SaveAs(Form("%s_%s_Cen%02d_%03d-%03d_eta.%s",sinputPrint.Data(),cAdd.Data(),ic,(Int_t)fStepLo,(Int_t)ptUp,printType.Data()));
327 c1->Modified();
328 c1->Update();
329 if(!gROOT->IsBatch())
330 if(getchar()=='q')return 1;
331 fStepLo += fStep;
332 fStepUp += fStep;
333 delete hProjAll;
334 delete hProjLead;
335 delete hProjLeadFull;
336 }
337
338 // Phi vs. p_T
339
340
341 TH2F *h2PhiPtLead[nJets+2] = {0,};
342 TH2F *h2PhiPtLeadFull[nJets+2] = {0,};
343 TH2F *h2PhiPtAll[nJets+2] = {0,};
344 HistsFromSingleJets(Form("fh2PhiPt%s_j%%d",cAdd.Data()),inputList,1,h2PhiPtLead);
345 HistsFromSingleJets(Form("fh2PhiPt%s_j%%d",cAddFull.Data()),inputList,1,h2PhiPtLeadFull);
346 HistsFromSingleJets(Form("fh2PhiPt%s_j%%d",cAddFull.Data()),inputList,nJets+1,h2PhiPtAll,1,nJets);
347
348 c1->SetLogz();
349 h2PhiPtLead[0]->Draw("colz");
350 c1->SaveAs(Form("%s_%s_Cen%02d_phiPtLead2D.%s",sinputPrint.Data(),cAdd.Data(),ic,printType.Data()));
351 c1->Modified();
352 c1->Update();
353
354 if(!gROOT->IsBatch())
355 if(getchar()=='q')return 1;
356
357 fStepLo = 0;
358 fStepUp = fStepLo+fStep;
359
360 while(fStepUp<fMax){
361 Int_t iBinLo = h2PhiPtLead[0]->GetYaxis()->FindBin(fStepLo);
362 Int_t iBinUp = h2PhiPtLead[0]->GetYaxis()->FindBin(fStepUp)-1;
363 TH1D *hProjLead = h2PhiPtLead[0]->ProjectionX(Form("hPhiPtLead_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
364 TH1D *hProjLeadFull = h2PhiPtLeadFull[0]->ProjectionX(Form("hPhiPtLeadFull_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
365 TH1D *hProjAll = h2PhiPtAll[nJets+1]->ProjectionX(Form("hPhiPtAll_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
366
367 ScaleH1(hProjLead,"#phi","1/N_{evt} dN_{jet}/d#phi",1./fNevents,true);
368 SetStyleH1(hProjLead,kBlue,0,kFullCircle);
369 ScaleH1(hProjLeadFull,"#phi","1/N_{evt} dN_{jet}/d#phi",1./fNevents,true);
370 SetStyleH1(hProjLeadFull,kRed,0,kFullCircle);
371 ScaleH1(hProjAll,"#phi","1/N_{evt} dN_{jet}/d#phi",1./fNevents,true);
372 SetStyleH1(hProjAll,kGray,1001,kFullCircle);
373
374 hProjAll->SetMinimum(0);
375
376 hProjAll->DrawCopy("hist");
377 hProjLead->DrawCopy("histsame");
378 hProjLeadFull->DrawCopy("histsame");
379
380 Float_t ptUp = fStepUp-h2PhiPtLead[0]->GetYaxis()->GetBinWidth(iBinUp);
381
382 Printf("%3.0f-%3.0f",fStepLo,ptUp);
383 TLatex *txt = 0;
384 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));
385 txt->SetTextSize(gStyle->GetTextSize());
386 txt->Draw();
387 c1->SaveAs(Form("%s_%s_Cen%02d_%03d-%03d_phi.%s",sinputPrint.Data(),cAdd.Data(),ic,(Int_t)fStepLo,(Int_t)ptUp,printType.Data()));
388 c1->Modified();
389 c1->Update();
390 if(!gROOT->IsBatch())
391 if(getchar()=='q')return 1;
392 fStepLo += fStep;
393 fStepUp += fStep;
394 delete hProjAll;
395 delete hProjLead;
396 delete hProjLeadFull;
397 }
398
399 // area vs. p_T
400
401
402 TH2F *h2AreaPtLead[nJets+2] = {0,};
403 TH2F *h2AreaPtLeadFull[nJets+2] = {0,};
404 TH2F *h2AreaPtAll[nJets+2] = {0,};
e2081c22 405 HistsFromSingleJets(Form("fh2AreaPt%s_j%%d",cAdd.Data()),inputList,1,h2AreaPtLead);
406 HistsFromSingleJets(Form("fh2AreaPt%s_j%%d",cAddFull.Data()),inputList,1,h2AreaPtLeadFull);
407 HistsFromSingleJets(Form("fh2AreaPt%s_j%%d",cAddFull.Data()),inputList,nJets+1,h2AreaPtAll,1,nJets);
408
409 c1->SetLogz();
410 h2AreaPtLead[0]->Draw("colz");
411 c1->Modified();
412 c1->Update();
413 c1->SaveAs(Form("%s_%s_Cen%02d_areaPtLead2D.%s",sinputPrint.Data(),cAdd.Data(),ic,printType.Data()));
414 if(!gROOT->IsBatch())
415 if(getchar()=='q')return 1;
416
417 fStepLo = 0;
418 fStepUp = fStepLo+fStep;
419
420 while(fStepUp<fMax){
421 Int_t iBinLo = h2AreaPtLead[0]->GetYaxis()->FindBin(fStepLo);
422 Int_t iBinUp = h2AreaPtLead[0]->GetYaxis()->FindBin(fStepUp)-1;
423 TH1D *hProjLead = h2AreaPtLead[0]->ProjectionX(Form("hAreaPtLead_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
424 TH1D *hProjLeadFull = h2AreaPtLeadFull[0]->ProjectionX(Form("hAreaPtLeadFull_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
425 TH1D *hProjAll = h2AreaPtAll[nJets+1]->ProjectionX(Form("hAreaPtAll_%d_%d_%d",iBinLo,iBinUp,ic),iBinLo,iBinUp,"e");
426
427 ScaleH1(hProjLead,"Area A","1/N_{evt} dN_{jet}/dA",1./fNevents,true);
428 SetStyleH1(hProjLead,kBlue,0,kFullCircle);
429 ScaleH1(hProjLeadFull,"Area A","1/N_{evt} dN_{jet}/dA",1./fNevents,true);
430 SetStyleH1(hProjLeadFull,kRed,0,kFullCircle);
431 ScaleH1(hProjAll,"Area A","1/N_{evt} dN_{jet}/dA",1./fNevents,true);
432 SetStyleH1(hProjAll,kGray,1001,kFullCircle);
433
434 hProjAll->SetMinimum(0);
435 hProjAll->DrawCopy("hist");
436 hProjLead->DrawCopy("histsame");
437 hProjLeadFull->DrawCopy("histsame");
438
439 Float_t ptUp = fStepUp-h2AreaPtLead[0]->GetYaxis()->GetBinWidth(iBinUp);
440
441 Printf("%3.0f-%3.0f",fStepLo,ptUp);
442 TLatex *txt = 0;
443 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));
444 txt->SetTextSize(gStyle->GetTextSize());
445 txt->Draw();
446 c1->SaveAs(Form("%s_%s_Cen%02d_%03d-%03d_area.%s",sinputPrint.Data(),cAdd.Data(),ic,(Int_t)fStepLo,(Int_t)ptUp,printType.Data()));
447 c1->Modified();
448 c1->Update();
449 if(!gROOT->IsBatch())
450 if(getchar()=='q')return 1;
451 fStepLo += fStep;
452 fStepUp += fStep;
453 delete hProjAll;
454 delete hProjLead;
455 delete hProjLeadFull;
456 }
457
458 // Area vs eta:
459 TH2F *h2AreaEta[nJets+1] = {0,};
460 HistsFromSingleJets(Form("fh2EtaArea%s_j%%d",cAdd.Data()),inputList,nJets,h2AreaEta);
461 ScaleH2(h2AreaEta[nJets],"area A","#eta","1/N_{evt} dN_{jets}/dAd#eta",1./fNevents,true);
462 c1->SetLogz();
463 h2AreaEta[nJets]->Draw("colz");
464 c1->Modified();
465 c1->Update();
466 c1->SaveAs(Form("%s_%s_Cen%02d_areaEta2D.%s",sinputPrint.Data(),cAdd.Data(),ic,printType.Data()));
467 if(!gROOT->IsBatch())
468 if(getchar()=='q')return 1;
469 }
470
471 return 0;
472}
473
5f5e3e77 474Int_t PlotSpectraPbPb2(UInt_t iPlotFlag,UInt_t iCenFlag,UInt_t iSpecFlag,Int_t iALICEType){
8dc333da 475
476 // Using now the THNSparse histos
477 //
478 const int kMaxFiles = 2;
479 TString sinputFile[kMaxFiles];
480 TString sinputDir[kMaxFiles];
8f206bf0 481 TString sinputJet[kMaxFiles];
8dc333da 482
8f206bf0 483 TCanvas *c1 = new TCanvas("c1","",800,600);
8dc333da 484 c1->SetLogy();
485 THnSparseF *fhnJetPtRec[kMaxFiles] = {0,};
486 THnSparseF *fhnTrackPtRec[kMaxFiles] = {0,};
487 THnSparseF *fhnEvent[kMaxFiles] = {0,};
8f206bf0 488 TH1F *fh1Centrality[kMaxFiles] = {0,};
8dc333da 489
490 const int nCen = 4;
491 const Float_t nColl[nCen] = {1502,742.5,250,46.6};
492 const Float_t fCentLo[nCen] = {0,10,30,50};
493 const Float_t fCentUp[nCen] = {10,30,50,80};
8f206bf0 494 const Int_t iColCen[nCen] = {kRed,kOrange+1,kBlue,kGray+1};
495
496 const Int_t iMarkerF[kMaxFiles] = {kFullCircle,kOpenCircle};// for putting different jet cents onto one plot
497
498 TH1D* hSpectrumTracks[kMaxFiles][nCen] = {0,};
499 TH1D* hSpectrumJets[kMaxFiles][nCen] = {0,};
500
501 bool scaleNevent = true;
8dc333da 502
503 Int_t iF = 0;
8f206bf0 504 sinputFile[iF] = "~/Dropbox/SharedJets/Christian/Files/PWG4_JetTasksOutput_Merge_LHC10h_AOD_110429a.root";
8dc333da 505 // sinputDir[iF] = "PWG4_spec2_clustersAOD_ANTIKT04_B2_Filter00128_Cut02000_Skip02_clustersAOD_ANTIKT04_B0_Filter00128_Cut02000_Skip02_0000000000_Class00";
506 sinputDir[iF] = "PWG4_spec2_clustersAOD_ANTIKT04_B2_Filter00128_Cut00150_Skip02__0000000000_Class00";
8f206bf0 507 sinputJet[iF] = "anti-k_{T}, R = 0.4 (150 MeV cut)";
8dc333da 508 fhnJetPtRec[iF] = (THnSparseF*)GetObjectFromFile("fhnJetPtRec",sinputFile[iF].Data(),sinputDir[iF].Data());
509 fhnTrackPtRec[iF] = (THnSparseF*)GetObjectFromFile("fhnTrackPtRec",sinputFile[iF].Data(),sinputDir[iF].Data());
8f206bf0 510 fhnEvent[iF] = (THnSparseF*) GetObjectFromFile("fhnEvent",sinputFile[iF].Data(),sinputDir[iF].Data());
511 fh1Centrality[iF] = (TH1F*)GetObjectFromFile("fh1Centrality",sinputFile[iF].Data(),sinputDir[iF].Data());
8dc333da 512
8f206bf0 513 Int_t iF = 1;
514 sinputFile[iF] = "~/Dropbox/SharedJets/Christian/Files/PWG4_JetTasksOutput_Merge_LHC10h_AOD_110429a.root";
515 sinputDir[iF] = "PWG4_spec2_clustersAOD_ANTIKT04_B2_Filter00128_Cut02000_Skip02_clustersAOD_ANTIKT04_B0_Filter00128_Cut02000_Skip02_0000000000_Class00";
516 sinputJet[iF] = "anti-k_{T}, R = 0.4 (2 GeV cut)";
517 fhnJetPtRec[iF] = (THnSparseF*)GetObjectFromFile("fhnJetPtRec",sinputFile[iF].Data(),sinputDir[iF].Data());
518 fhnTrackPtRec[iF] = (THnSparseF*)GetObjectFromFile("fhnTrackPtRec",sinputFile[iF].Data(),sinputDir[iF].Data());
519 fh1Centrality[iF] = (TH1F*)GetObjectFromFile("fh1Centrality",sinputFile[iF].Data(),sinputDir[iF].Data());
8dc333da 520
521 fhnEvent[iF] = (THnSparseF*) GetObjectFromFile("fhnEvent",sinputFile[iF].Data(),sinputDir[iF].Data());
522
8f206bf0 523 TLatex *txtHead = new TLatex();
524 txtHead->SetNDC();
525 txtHead->SetTextAlign(23);
526
527
528 Float_t yMin = 1;
8dc333da 529
8f206bf0 530 for(int iF = 0;iF<kMaxFiles;iF++){
531 Bool_t bFirst1 = true;
532 for(int ic = 0;ic <nCen;ic++){
533 // project out the jets
534 Int_t ibLo = fhnJetPtRec[iF]->GetAxis(2)->FindBin(fCentLo[ic]+0.001);
535 Int_t ibUp = fhnJetPtRec[iF]->GetAxis(2)->FindBin(fCentUp[ic]-0.001);
8dc333da 536
8f206bf0 537 Float_t nEvents = fh1Centrality[iF]->Integral(fh1Centrality[iF]->FindBin(fCentLo[ic]+0.001),
538 fh1Centrality[iF]->FindBin(fCentUp[ic]-0.001));
539
540 fhnJetPtRec[iF]->GetAxis(2)->SetRange(ibLo,ibUp);
541 fhnJetPtRec[iF]->GetAxis(0)->SetRange(3,3); // take all jets
542
543 hSpectrumJets[iF][ic] = fhnJetPtRec[iF]->Projection(1,"E");
544 hSpectrumJets[iF][ic]->SetName(Form("hSpectrumJets_%d_C%d",iF,ic));
545
546 SetHistoAttributes(hSpectrumJets[iF][ic],kFullCircle,iColCen[ic]);
547 hSpectrumJets[iF][ic]->Rebin(2);
5f5e3e77 548 hSpectrumJets[iF][ic]->SetAxisRange(0,100);
8f206bf0 549 hSpectrumJets[iF][ic]->SetMinimum(0.1);
550 hSpectrumJets[iF][ic]->SetXTitle("p_{T} (GeV/c)");
551
552
553 for(int ib = 1;ib <=hSpectrumJets[iF][ic]->GetNbinsX();ib++){
554 Float_t width = hSpectrumJets[iF][ic]->GetBinWidth(ib);
555 hSpectrumJets[iF][ic]->SetBinContent(ib,hSpectrumJets[iF][ic]->GetBinContent(ib)/width);
556 hSpectrumJets[iF][ic]->SetBinError(ib,hSpectrumJets[iF][ic]->GetBinError(ib)/width);
557 }
558 if(scaleNevent&&nEvents>0){
559 Printf("%d %E events",ic,nEvents);
560 hSpectrumJets[iF][ic]->Scale(1./nEvents);
561 hSpectrumJets[iF][ic]->SetYTitle("1/N_{evt} dN/dp_{T}d#eta");
562 if(yMin>1./nEvents&&(iCenFlag&1<<ic))yMin=1./nEvents;
563 }
564 else{
565 hSpectrumJets[iF][ic]->SetYTitle("dN/dp_{T}d#eta");
566 }
567
568 // project out the tracks
569 fhnTrackPtRec[iF]->GetAxis(2)->SetRange(ic+1,ic+1);
570 fhnTrackPtRec[iF]->GetAxis(0)->SetRange(2,2); // take all tracks
571
572 hSpectrumTracks[iF][ic] = fhnTrackPtRec[iF]->Projection(1,"E");
573 hSpectrumTracks[iF][ic]->SetName(Form("hSpectrumTracks_%d_C%d",iF,ic));
574
575 SetHistoAttributes(hSpectrumTracks[iF][ic],kFullCircle,iColCen[ic]);
576 hSpectrumTracks[iF][ic]->SetAxisRange(0,300);
577 hSpectrumTracks[iF][ic]->SetMinimum(0.1);
578 hSpectrumTracks[iF][ic]->SetXTitle("p_{T} (GeV/c)");
579
580
581 for(int ib = 1;ib <=hSpectrumTracks[iF][ic]->GetNbinsX();ib++){
582 Float_t width = hSpectrumTracks[iF][ic]->GetBinWidth(ib);
583 hSpectrumTracks[iF][ic]->SetBinContent(ib,hSpectrumTracks[iF][ic]->GetBinContent(ib)/width/1.8);
584 hSpectrumTracks[iF][ic]->SetBinError(ib,hSpectrumTracks[iF][ic]->GetBinError(ib)/width);
585 }
586 if(scaleNevent&&nEvents>0){
587 Printf("%d %E events",ic,nEvents);
588 hSpectrumTracks[iF][ic]->Scale(1./nEvents);
589 hSpectrumTracks[iF][ic]->SetYTitle("1/N_{evt} dN/dp_{T}d#eta");
590 }
591 else{
592 hSpectrumTracks[iF][ic]->SetYTitle("dN/dp_{T}d#eta");
593 }
594 }
595 }
596
597 // the drawin part
598
599 // draw the spectra for all cents
600 if(iPlotFlag&1<<0){
601
602 for(int iF = 0;iF<kMaxFiles;iF++){
603 Bool_t bFirst1 = true;
604
5f5e3e77 605 TLegend *leg1 = new TLegend(0.58,0.55,0.85,0.85);
8f206bf0 606 leg1->SetFillColor(0);
607 leg1->SetTextFont(gStyle->GetTextFont());
608 leg1->SetBorderSize(0);
609 leg1->SetTextAlign(12);
5f5e3e77 610 leg1->SetTextSize(0.028);
8f206bf0 611
5f5e3e77 612 if(!iSpecFlag&1){
613 leg1->SetHeader(sinputJet[iF].Data());
614 }
8f206bf0 615 for(int ic = 0;ic <nCen;ic++){
616 if(!(iCenFlag&(1<<ic)))continue;
617 c1->cd();
618 if(bFirst1){
619 hSpectrumJets[iF][ic]->SetMinimum(yMin/10.);
620 hSpectrumJets[iF][ic]->DrawCopy("p");
621 bFirst1 = false;
622 }
623 else {
624 hSpectrumJets[iF][ic]->DrawCopy("psame");
625 }
5f5e3e77 626 /*
8f206bf0 627 leg1->AddEntry((TObject*)0,Form("%2.0f - %2.0f%%",fCentLo[ic],fCentUp[ic]),"");
628 leg1->AddEntry(hSpectrumJets[iF][ic],Form("jets: %s",sinputJet[iF].Data()), "P");
5f5e3e77 629 */
8f206bf0 630 if(iSpecFlag&1){
631 hSpectrumTracks[iF][ic]->DrawCopy("histsame");
5f5e3e77 632 leg1->AddEntry((TObject*)0,Form("%2.0f - %2.0f%%",fCentLo[ic],fCentUp[ic]),"");
633 leg1->AddEntry(hSpectrumJets[iF][ic],Form("jets: %s",sinputJet[iF].Data()), "P");
8f206bf0 634 leg1->AddEntry(hSpectrumTracks[iF][ic],Form("tracks"), "L");
635 }
5f5e3e77 636 else{
637 leg1->AddEntry(hSpectrumJets[iF][ic],Form("%2.0f - %2.0f%%",fCentLo[ic],fCentUp[ic]), "P");
638 }
8f206bf0 639 }
640
641
642
8f206bf0 643
1354c5f3 644 txtHead->DrawLatex(0.5,0.99,"LHC2010 Pb-Pb #sqrt{s_{NN}} = 2.76 TeV");
8f206bf0 645 leg1->Draw();
5f5e3e77 646
647
648 DrawALICE(c1,iALICEType,0.3,0.35);
649
8f206bf0 650 c1->Update();
651 picName = Form("%sCentSpectraPbPb_Jets%s_%d_%d.%s",picPrefix.Data(),(iSpecFlag&1?"Tracks":""),iCenFlag,iF,picSuffix.Data());
652 c1->SaveAs(picName.Data());
5f5e3e77 653 if(!gROOT->IsBatch()){
654 if(getchar()=='q')return 1;
655 }
8dc333da 656 }
8f206bf0 657 }
8dc333da 658
8f206bf0 659 // draw the spectra for each cent togeter
660 if(iPlotFlag&1<<1){
661
662
663 for(int ic = 0;ic <nCen;ic++){
664
665 if(!(iCenFlag&(1<<ic)))continue;
666
667
668 Bool_t bFirst1 = true;
669
670 TLegend *leg1 = new TLegend(0.45,0.55,0.8,0.85);
671 leg1->SetFillColor(0);
672 leg1->SetTextFont(gStyle->GetTextFont());
673 leg1->SetBorderSize(0);
674 leg1->SetTextAlign(12);
675
676 for(int iF = 0;iF<kMaxFiles;iF++){
677 c1->cd();
678 hSpectrumJets[iF][ic]->SetMarkerStyle(iMarkerF[iF]);
679 if(bFirst1){
680 leg1->AddEntry((TObject*)0,Form("%2.0f - %2.0f%%",fCentLo[ic],fCentUp[ic]),"");
681 }
682 if(bFirst1){
683 hSpectrumJets[iF][ic]->SetMinimum(yMin/10.);
684 hSpectrumJets[iF][ic]->DrawCopy("p");
685 if(iSpecFlag&1){
686 hSpectrumTracks[iF][ic]->DrawCopy("histsame");
687 leg1->AddEntry(hSpectrumTracks[iF][ic],Form("tracks"), "L");
688 }
689 bFirst1 = false;
690 }
691 else {
692 hSpectrumJets[iF][ic]->DrawCopy("psame");
693 }
694 leg1->AddEntry(hSpectrumJets[iF][ic],Form("jets: %s",sinputJet[iF].Data()), "P");
695 }
1354c5f3 696 txtHead->DrawLatex(0.5,0.99,"LHC2010 Pb-Pb #sqrt{s_{NN}} = 2.76 TeV");
8f206bf0 697 leg1->Draw();
698
5f5e3e77 699 DrawALICE(c1,iALICEType,0.26,0.35);
8f206bf0 700 c1->Update();
701 picName = Form("%sCmpSpectraPbPb_Jets%s_%d.%s",picPrefix.Data(),(iSpecFlag&1?"Tracks":""),ic,picSuffix.Data());
702 c1->SaveAs(picName.Data());
5f5e3e77 703 if(!gROOT->IsBatch()){
704 if(getchar()=='q')return 1;
705 }
8f206bf0 706 }
8dc333da 707 }
8f206bf0 708
709
8dc333da 710 return 0;
711}
712
5f5e3e77 713Int_t PlotSubtract(Int_t iALICEType){
8dc333da 714
8f206bf0 715
716
717 gStyle->SetTitleYSize(0.7*gStyle->GetTitleYSize());
718 gStyle->SetTitleXSize(0.7*gStyle->GetTitleXSize());
719 gStyle->SetPadRightMargin(1.1*gStyle->GetPadRightMargin());
720 gStyle->SetPadLeftMargin(0.9*gStyle->GetPadLeftMargin());
721 gStyle->SetPadBottomMargin(0.7*gStyle->GetPadBottomMargin());
722 gROOT->ForceStyle();
723
724
8dc333da 725 // Using now the THNSparse histos
726 //
727
8f206bf0 728 TCanvas *c1 = new TCanvas("c11","",700,600); // create with good aspect ratios...
8dc333da 729 c1->SetLogz();
730
8f206bf0 731
8dc333da 732 const int kMaxFiles = 2;
733 TString sinputFile[kMaxFiles];
734 TString sinputDir[kMaxFiles];
735 TString sLegend[kMaxFiles];
736
737
738
739
740 // Parameter
741 Float_t fMaxMult[kMaxFiles] = {0,};
742 Float_t fMaxRho[kMaxFiles] = {0,};
8f206bf0 743 Float_t fMaxSigma[kMaxFiles] = {0,};
8dc333da 744
745
746 TH2F* fh2CentvsRho[kMaxFiles]= {0,};
747 TH2F* fh2MultvsRho[kMaxFiles]= {0,};
8f206bf0 748 TH2F* fh2CentvsSigma[kMaxFiles]= {0,};
749 TH2F* fh2MultvsSigma[kMaxFiles]= {0,};
8dc333da 750
751
752
753 const Int_t nSub = 3;
754 THnSparseF* hnDPtAreaCentMult[kMaxFiles][nSub] = {0,};
755
756 Int_t iF = 0;
8f206bf0 757 sinputFile[iF] = "~/Dropbox/SharedJets/Bastian/Files/PWG4_JetTasksOutput.root";
758 sinputDir[iF] = "pwg4JetSubtract_B2";
5f5e3e77 759 sLegend[iF] = "FastJet k_{T} R = 0.4 (cut 150 MeV)";
8dc333da 760 fMaxMult[iF] = 3500;
761 fMaxRho[iF] = 250;
8f206bf0 762 fMaxSigma[iF] = 40;
8dc333da 763 iF++;
764
8f206bf0 765 sinputFile[iF] = "~/Dropbox/SharedJets/Bastian/Files/PWG4_JetTasksOutput.root";
766 sinputDir[iF] = "pwg4JetSubtract_B2_Cut2000";
8dc333da 767 sLegend[iF] = "B2 (2000 MeV)";
768 fMaxMult[iF] = 300;
769 fMaxRho[iF] = 50;
8f206bf0 770 fMaxSigma[iF] = 10;
8dc333da 771
772
773
774 for(iF = 0;iF<kMaxFiles;iF++){
775 if(sinputFile[iF].Length()==0)continue;
776 // fetch all the histos
8f206bf0 777 fh2CentvsRho[iF]= (TH2F*) GetObjectFromFile("fh2CentvsRho",sinputFile[iF].Data(),sinputDir[iF].Data(),"PWG4_","pwg4");
778 fh2MultvsRho[iF]= (TH2F*) GetObjectFromFile("fh2MultvsRho",sinputFile[iF].Data(),sinputDir[iF].Data(),"PWG4_","pwg4");
779 fh2CentvsSigma[iF]= (TH2F*) GetObjectFromFile("fh2CentvsSigma",sinputFile[iF].Data(),sinputDir[iF].Data(),"PWG4_","pwg4");
780 fh2MultvsSigma[iF]= (TH2F*) GetObjectFromFile("fh2MultvsSigma",sinputFile[iF].Data(),sinputDir[iF].Data(),"PWG4_","pwg4");
781
782
783 if(!fh2CentvsRho[iF])continue;
8dc333da 784
785 for(int iSub = 0;iSub<nSub;iSub++){
786 hnDPtAreaCentMult[iF][iSub] = (THnSparseF*) GetObjectFromFile(Form("hnDPtAreaCentMult_%d",iSub),sinputFile[iF].Data(),sinputDir[iF].Data(),"pwg");
787 }
8f206bf0 788
789 SetHistAxisStyle(fh2CentvsRho[iF],1.5,1.2,1.1,
790 0,1.1,1.1);
791 SetHistAxisStyle(fh2MultvsRho[iF],1.5,1.2,1.1,
792 0,1.1,1.1);
793
794 SetHistAxisStyle(fh2CentvsSigma[iF],1.5,1.2,1.1,
795 0,1.1,1.1);
796 SetHistAxisStyle(fh2MultvsSigma[iF],1.5,1.2,1.1,
797 0,1.1,1.1);
798
799
800
8dc333da 801 TLatex *txt = new TLatex();
802 txt->SetNDC();
8f206bf0 803 txt->SetTextSize(gStyle->GetTextSize()*0.8);
5f5e3e77 804 txt->SetTextAlign(22);
8f206bf0 805
806 TLatex *txtHead = new TLatex();
807 txtHead->SetNDC();
808 txtHead->SetTextAlign(23);
809
1354c5f3 810 TProfile *hProf = fh2CentvsRho[iF]->ProfileY(Form("%s_profile_%d",fh2CentvsRho[iF]->GetName(),iF,1,-1,"s"));
811
812 fh2CentvsRho[iF]->SetAxisRange(0,400,"Y");
813 TGraphErrors *grRho = new TGraphErrors(fh2CentvsRho[iF]->GetNbinsX());
814 grRho->SetLineColor(kBlack);
815 for(int ib=1;ib<=fh2CentvsRho[iF]->GetNbinsX();ib++){
816 TH1D *hProj = fh2CentvsRho[iF]->ProjectionY("hProj",ib,ib);
817 hProj->SetBinContent(1,0);
818 hProj->SetBinError(1,0);
819 Float_t mean = hProj->GetMean();
820 Float_t cent = fh2CentvsRho[iF]->GetXaxis()->GetBinCenter(ib);
821 Float_t rms = hProj->GetRMS();
822 Float_t nEntries = hProj->Integral();
823 grRho->SetPoint(ib-1,cent,mean);
824 // if(nEntries)grRho->SetPointError(ib-1,0,rms/TMath::Sqrt(nEntries));
825 hProj->Delete();
826 }
827
8f206bf0 828
8dc333da 829
830 fh2CentvsRho[iF]->SetAxisRange(0.,fMaxRho[iF],"Y");
8f206bf0 831 fh2CentvsRho[iF]->SetYTitle("#rho (GeV/c)");
832 fh2CentvsRho[iF]->SetXTitle("V0 centrality (%)");
833 fh2CentvsRho[iF]->GetXaxis()->SetRangeUser(0,80);
8dc333da 834 fh2CentvsRho[iF]->Draw("colz");
1354c5f3 835
836 grRho->Draw("sameC");
837
5f5e3e77 838 txt->DrawLatex(0.5,0.8,sLegend[iF].Data());
8f206bf0 839 txtHead->DrawLatex(0.5,0.99,"LHC2010 PbPb #sqrt{s_{NN}} = 2.76 TeV");
840
5f5e3e77 841 DrawALICE(c1,iALICEType,0.6,0.6);
8dc333da 842 c1->Update();
843 picName = Form("%sCentVsRho_%d.%s",picPrefix.Data(),iF,picSuffix.Data());
844 c1->SaveAs(picName.Data());
5f5e3e77 845 if(!gROOT->IsBatch()){
846 if(getchar()=='q')return 1;
847 }
1354c5f3 848
8dc333da 849 fh2MultvsRho[iF]->SetAxisRange(0.,fMaxRho[iF],"Y");
850 fh2MultvsRho[iF]->SetAxisRange(0.,fMaxMult[iF]);
8f206bf0 851 fh2MultvsRho[iF]->SetXTitle("N^{raw}_{input}");
852 fh2MultvsRho[iF]->SetYTitle("#rho (GeV/c)");
853
1354c5f3 854 fh2CentvsRho[iF]->SetAxisRange(0,400,"Y");
855 Int_t ibMaxM = fh2MultvsRho[iF]->GetXaxis()->FindBin(2600);
856 TGraphErrors *grRhoM = new TGraphErrors(ibMaxM-6);
857 grRhoM->SetLineColor(kBlack);
858 for(int ib=1;ib<ibMaxM;ib++){
859 if(ib<4)continue;
860 TH1D *hProj = fh2MultvsRho[iF]->ProjectionY("hProj",ib,ib);
861 hProj->SetBinContent(1,0);
862 hProj->SetBinError(1,0);
863 Float_t mean = hProj->GetMean();
864 Float_t cent = fh2MultvsRho[iF]->GetXaxis()->GetBinCenter(ib);
865 Float_t rms = hProj->GetRMS();
866 Float_t nEntries = hProj->Integral();
867 if(nEntries>0){
868 grRhoM->SetPoint(ib-6,cent,mean);
869 // grRhoM->SetPointError(ib-6,0,rms/TMath::Sqrt(nEntries));
870 }
871 hProj->Delete();
872 }
873 grRhoM->Print();
874
8dc333da 875 fh2MultvsRho[iF]->Draw("colz");
1354c5f3 876 grRhoM->Draw("Csame");
5f5e3e77 877 txt->DrawLatex(0.5,0.8,sLegend[iF].Data());
8f206bf0 878 txtHead->DrawLatex(0.5,0.99,"LHC2010 PbPb #sqrt{s_{NN}} = 2.76 TeV");
5f5e3e77 879 DrawALICE(c1,iALICEType,0.4,0.7);
8dc333da 880 c1->Update();
881 picName = Form("%sMultVsRho_%d.%s",picPrefix.Data(),iF,picSuffix.Data());
882 c1->SaveAs(picName.Data());
5f5e3e77 883 if(!gROOT->IsBatch()){
884 if(getchar()=='q')return 1;
885 }
8f206bf0 886
887 fh2MultvsSigma[iF]->SetAxisRange(0.,fMaxMult[iF]);
888 fh2MultvsSigma[iF]->SetAxisRange(0.,fMaxSigma[iF],"Y");
889 fh2MultvsSigma[iF]->SetXTitle("N^{raw}_{input}");
890 fh2MultvsSigma[iF]->SetYTitle("#sigma (GeV/c)");
891
892 fh2MultvsSigma[iF]->Draw("colz");
893 txt->DrawLatex(0.3,0.7,sLegend[iF].Data());
894 txtHead->DrawLatex(0.5,0.99,"LHC2010 PbPb #sqrt{s_{NN}} = 2.76 TeV");
5f5e3e77 895 DrawALICE(c1,iALICEType);
8f206bf0 896 c1->Update();
897 picName = Form("%sMultVsSigma_%d.%s",picPrefix.Data(),iF,picSuffix.Data());
898 c1->SaveAs(picName.Data());
5f5e3e77 899 if(!gROOT->IsBatch()){
900 if(getchar()=='q')return 1;
901 }
8f206bf0 902
903 fh2CentvsSigma[iF]->SetAxisRange(0.,fMaxSigma[iF],"Y");
904 fh2CentvsSigma[iF]->SetYTitle("#sigma (GeV/c)");
905 fh2CentvsSigma[iF]->SetXTitle("V0 centrality (%)");
906 fh2CentvsSigma[iF]->GetXaxis()->SetRangeUser(0,80);
907 fh2CentvsSigma[iF]->Draw("colz");
908
909
910 fh2CentvsSigma[iF]->Draw("colz");
911 txt->DrawLatex(0.3,0.7,sLegend[iF].Data());
912 txtHead->DrawLatex(0.5,0.99,"LHC2010 PbPb #sqrt{s_{NN}} = 2.76 TeV");
5f5e3e77 913 DrawALICE(c1,iALICEType);
8f206bf0 914 c1->Update();
915 picName = Form("%sCentVsSigma_%d.%s",picPrefix.Data(),iF,picSuffix.Data());
916 c1->SaveAs(picName.Data());
5f5e3e77 917 if(!gROOT->IsBatch()){
918 if(getchar()=='q')return 1;
919 }
8f206bf0 920
921
8dc333da 922 }
923
924 return 0;
925 Bool_t bFirst1 = true;
926 iF = 0;
927 for(int ic = 0;ic <10;ic++){
928 Int_t iSub = 2;
929 hnDPtAreaCentMult[iF][iSub]->GetAxis(2)->SetRange(ic+1,ic+1);
930 TH2D *hDptArea = hnDPtAreaCentMult[iF][iSub]->Projection(0,1,"E");
931 hDptArea->SetName(Form("hDptArea_C%d_%d",ic,iSub));
932 hDptArea->Draw("colz");
933 c1->Update();
934 if(getchar()=='q')return 1;
935 }
936
937
938 return 0;
939
940}
941
942
e2081c22 943Int_t PlotSpectraPbPb(){
944
945 // PLot the simple 1D histos from the spectrum task
946 //
947 const Int_t iRebin = 4;
948 const int kMaxFiles = 5;
949
950 Double_t yMin = 0.01;
951 Double_t yMax = 1E+07;
952
953 TString sinputFile[kMaxFiles];
954 TFile* inputFile[kMaxFiles] = {kMaxFiles*0};
955 TString sinputDir[kMaxFiles];
956 TDirectory *inputDir[kMaxFiles] = {kMaxFiles*0};
957 TString sinputList[kMaxFiles];
958 TString sCen[kMaxFiles];
959 TList* inputList[kMaxFiles] = {kMaxFiles*0};
960 // TString sinputLegend[kMaxFiles];
961 TString sinputJet[kMaxFiles];
962 Int_t kColor[kMaxFiles];
963 Float_t nColl[kMaxFiles];
964
965 bool bLogLog = false;
966 // const Int_t kRef = 1; // anti -kT
967
968 //
969 TH1F* hJets[kMaxFiles];
970
971 const int nJets = 2;
972
973 Int_t iJF = 0;
974 sinputFile[iJF] = "~/alice/data/analysis/train/LHC10h_110116/PWG4_JetTasksOutput.root";
975 // sinputFile[iJF] = "~/alice/data/analysis/train_maf/output_110210.root";
976 sinputDir[iJF] = "PWG4_spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class00";
977 sinputList[iJF] = "pwg4spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class00";
978 sinputJet[iJF] = "anti-k_{T} R = 0.4";
979 sCen[iJF] = " 0-80%";
980 kColor[iJF] = kBlack;
981 nColl[iJF] = 453.3;
982 iJF++;
983
984 // sinputFile[iJF] = "../output/PWG4_JetTasksOutput_Merge_051205b_Subtract3_UA1.root";
985 sinputFile[iJF] = "~/alice/data/analysis/train/LHC10h_110116/PWG4_JetTasksOutput.root";
986 // sinputFile[iJF] = "~/alice/data/analysis/train_maf/output_110210.root";
987 sinputDir[iJF] = "PWG4_spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class01";
988 sinputList[iJF] = "pwg4spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class01";
989 sinputJet[iJF] = " 0-10% Anti-k_{T} R = 0.4";
990 sCen[iJF] = " 0-10%";
991 kColor[iJF] = kRed;
992 nColl[iJF] = 1502;
993 iJF++;
994
995 // sinputFile[iJF] = "../output/PWG4_JetTasksOutput_Merge_051205b_Subtract3_UA1.root";
996 sinputFile[iJF] = "~/alice/data/analysis/train/LHC10h_110116/PWG4_JetTasksOutput.root";
997 // sinputFile[iJF] = "~/alice/data/analysis/train_maf/output_110210.root";
998 sinputDir[iJF] = "PWG4_spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class02";
999 sinputList[iJF] = "pwg4spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class02";
1000 sinputJet[iJF] = "10-30% Anti-k_{T} R = 0.4";
1001 sCen[iJF] = "10-30%";
1002 kColor[iJF] = kRed-4;
1003 nColl[iJF] = 742.5;
1004 iJF++;
1005
1006 // sinputFile[iJF] = "../output/PWG4_JetTasksOutput_Merge_051205b_Subtract3_UA1.root";
1007 sinputFile[iJF] = "~/alice/data/analysis/train/LHC10h_110116/PWG4_JetTasksOutput.root";
1008 // sinputFile[iJF] = "~/alice/data/analysis/train_maf/output_110210.root";
1009 sinputDir[iJF] = "PWG4_spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class03";
1010 sinputList[iJF] = "pwg4spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class03";
1011 sinputJet[iJF] = "30-50% Anti-k_{T} R = 0.4";
1012 sCen[iJF] = "30-50%";
1013 kColor[iJF] = kRed-7;
1014 nColl[iJF] = 250;
1015 iJF++;
1016
1017 // sinputFile[iJF] = "../output/PWG4_JetTasksOutput_Merge_051205b_Subtract3_UA1.root";
1018 sinputFile[iJF] = "~/alice/data/analysis/train/LHC10h_110116/PWG4_JetTasksOutput.root";
1019 // sinputFile[iJF] = "~/alice/data/analysis/train_maf/output_110210.root";
1020 sinputDir[iJF] = "PWG4_spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class04";
1021 sinputList[iJF] = "pwg4spec2_clustersAOD_ANTIKT04_B3_Filter00256_Cut00150_Skip00_clustersAOD_ANTIKT04_B1_Filter00256_Cut00150_Skip00_0000000000_Class04";
1022 sinputJet[iJF] = "Anti-k_{T} R = 0.4";
1023 sCen[iJF] = "50-80%";
1024 kColor[iJF] = kRed-9;
1025 nColl[iJF] = 46.6;
1026 iJF++;
1027
1028
1029 for(int iF = 0; iF < kMaxFiles;iF++){
1030 if(sinputFile[iF].Length()==0)continue;
1031 if(gROOT->GetListOfFiles()->FindObject(sinputFile[iF].Data())){
1032 inputFile[iF] = (TFile*)gROOT->GetListOfFiles()->FindObject(sinputFile[iF].Data());
1033 }
1034 else {
1035 inputFile[iF] = TFile::Open(sinputFile[iF].Data());
1036 }
1037 inputDir[iF] = (TDirectory*)inputFile[iF]->Get(sinputDir[iF].Data());
1038 if(!inputDir[iF]){
1039 Printf("Dir not found %s",sinputDir[iF].Data());
1040 continue;
1041 }
1042 inputList[iF] = (TList*)inputDir[iF]->Get(sinputList[iF].Data());
1043 if(!inputList[iF]){
1044 Printf("List not found %s",sinputList[iF].Data());
1045 continue;
1046 }
1047 }
1048
1049
1050
1051
1052 TCanvas *c1 = new TCanvas("c1","spectra",20,20,800,600);
1053 TCanvas *c2 = new TCanvas("c2","ratio",820,20,800,600);
1054
1055 TLegend *leg1 = new TLegend(0.45,0.65,0.7,0.85);
1056 leg1->SetHeader(Form("%s |#eta| < 0.4",sinputJet[0].Data()));
1057 leg1->SetFillColor(0);
1058 leg1->SetTextFont(gStyle->GetTextFont());
1059 leg1->SetBorderSize(0);
1060
1061 char *cMaskRec = "fh1PtRecIn_j%d";
1062
1063 c1->SetLogy();
1064 if(bLogLog)c1->SetLogx();
1065
1066 for(int iF = 0; iF < kMaxFiles;iF++){
1067 if(sinputFile[iF].Length()==0)continue;
1068 // Draw the jet spectra
1069 if(!inputList[iF])continue;
1070 TH1F *hRec[nJets+1];
1071 HistsFromSingleJets(cMaskRec,inputList[iF],nJets,hRec,iRebin);
1072 /// TH1F *hRec[nJets] = (TH1F*)inputList[iF]->FindObject("fh1PtJetsRecIn"); hRec[nJets]->Rebin(iRebin);
1073
1074 hRec[nJets]->SetMarkerStyle(kFullCircle);
1075 hRec[nJets]->SetXTitle("p_{T} (GeV/c)");
1076 hRec[nJets]->SetMarkerColor(kColor[iF]);
1077 hRec[nJets]->Scale(1./1.6); // deta
1078 hRec[nJets]->Scale(1./hRec[nJets]->GetBinWidth(1));
1079
1080 hRec[nJets]->SetYTitle("dN/dp_{T}dy");
1081
1082 if(true){
1083
1084 ScaleNevent(hRec[nJets], inputFile[iF],AliAnalysisTaskJetServices::kSelected,0,iF); // <----- careful assumes iCen == iF
1085 hRec[nJets]->SetYTitle("1/N_{evt} dN/dp_{T}");
1086 // Printf("%d Ncall %f",iF,nColl[iF]);
1087 // hRec[nJets]->Scale(1./nColl[iF]);
1088 // hRec[nJets]->SetYTitle("1/(N_{coll}*(N_{evt}) dN/dp_{T}dy");
1089 yMin = 1E-8;
1090 yMax = 1;
1091 }
1092 if(sCen[iF].Length()){
1093 leg1->AddEntry(hRec[nJets],Form("%s",sCen[iF].Data()),"P");
1094 }
1095
1096 hJets[iF] = hRec[nJets];
1097 c1->cd();
1098
1099
1100 if(bLogLog)hRec[nJets]->SetAxisRange(5,100);
1101 else hRec[nJets]->SetAxisRange(1,100);
1102
1103 if(iF==0){
1104 if(yMax>0) hRec[nJets]->SetMaximum(yMax);
1105 if(yMin>0) hRec[nJets]->SetMinimum(yMin);
1106 hRec[nJets]->DrawCopy("P");
1107 }
1108 else{
1109 if(sinputJet[iF].Length())hRec[nJets]->DrawCopy("Psame");
1110 }
1111 c1->Update();
1112 if(getchar()=='q')return 1;
1113 }
1114
1115
1116 c1->Update();
8dc333da 1117 picName = "jetspectrumPbPb";
e2081c22 1118 if(bLogLog)picName += "LogLog";
1119 leg1->Draw("same");
1120 TLatex *txt = 0;
1354c5f3 1121 txt = new TLatex(5,0.1,"LHC2010 Pb-Pb #sqrt{s_{NN}} = 2.76 TeV");
e2081c22 1122 txt->SetTextSize(gStyle->GetTextSize()*0.8);
1123 txt->Draw();
1124 ALICEWorkInProgress(c1,"02/15/2011");
1125 c1->SaveAs(Form(cPrintMask,picName.Data()));
1126 if(getchar()=='q')return 1;
1127
1128
1129 c2->cd();
1130 leg1->Draw("same");
1131
1132
1133
1134
1135 picName = "jetspectrumLabels";
1136 c2->SaveAs(Form(cPrintMask,picName.Data()));
1137 if(getchar()=='q')return 1;
1138
1139 // scale with nColl
1140 for(int iF = 0; iF < kMaxFiles;iF++){
1141 if(sinputFile[iF].Length()==0)continue;
1142 if(!hJets[iF])continue;
1143 Printf("%d Ncoll %f",iF,nColl[iF]);
1144 hJets[iF]->Scale(1./nColl[iF]);
1145 hJets[iF]->SetYTitle("1/(N_{coll}*(N_{evt}) dN/dp_{T}dy");
1146 }
1147
1148 const Int_t nPSmear = 9;
1149 TH1F *hPythia[nPSmear];
1150 // Fetch the smeared pythia spectra
1151
1152 TString pName;
1153 c2->SetLogy();
1154 TH1F* hRatio[kMaxFiles];
1155
1156 for(int i = 0;i < nPSmear;i++){
1157 if(i==nPSmear-1){
1158 pName = "h1JetPtCh";
1159 hPythia[i] = GetPythia8Spectrum(pName.Data(),"~/alice/sim/pythia8/examples/output/LHCJets_2.75TeV_allpt_R04.root",1.0,iRebin);
1160 }
1161 else {
1162 pName = Form("h1JetPtCh_C%02d_j%%d",i);
1163 hPythia[i] = GetPythia8Spectrum(pName.Data(),"~/alice/sim/pythia8/examples/output/LHCJets_2.75TeV_allpt_R04.root",1.0,iRebin,nJets);
1164 }
1165 hPythia[i]->SetMarkerStyle(kOpenCircle);
1166 hPythia[i]->SetMarkerColor(kBlack);
1167 hPythia[i]->SetYTitle("1/(N_{coll}*N_{evt}) dN/dp_{T}dy");
1168 hPythia[i]->SetAxisRange(0,100);
1169 c1->cd();
1170 hPythia[i]->Draw("CP");
1171 hPythia[i]->SetMaximum(1E-02);
1172 hPythia[i]->SetMinimum(1E-09);
1173 for(int iF = 0; iF < kMaxFiles;iF++){
1174 if(sinputFile[iF].Length()==0)continue;
1175 if(!hJets[iF])continue;
1176
1177 c1->cd();
1178 hJets[iF]->Draw("psame");
1179
1180 c2->cd();
1181 hRatio[iF] = (TH1F*)hJets[iF]->Clone(Form("hRatio%d",iF));
1182 hRatio[iF]->SetMaximum(100);
1183 hRatio[iF]->SetMinimum(0.01);
1184 hRatio[iF]->Divide(hPythia[i]);
1185 hRatio[iF]->SetYTitle("Ratio");
1186 if(iF==0) hRatio[iF]->DrawCopy("p");
1187 else hRatio[iF]->DrawCopy("psame");
1188 c1->Update();
1189 c2->Update();
1190 }
1191
1192 if(getchar()=='q')return 1;
1193 picName = Form("jetspectrumRatioPbPb_Smear%d",i);
1194 c2->SaveAs(Form(cPrintMask,picName.Data()));
1195 picName = Form("jetspectrumPbPb_Smear%d",i);
1196 c1->SaveAs(Form(cPrintMask,picName.Data()));
1197 if(getchar()=='q')return 1;
1198
1199 }
1200 c1->cd();
1201 hPythia[nPSmear-1]->DrawCopy("P");
1202 for(int i = 0;i < nPSmear-1;i++){
1203 c1->cd();
1204 hPythia[i]->SetMarkerStyle(kFullCircle);
1205 hPythia[i]->SetMarkerColor(kRed-9+i);
1206 hPythia[i]->DrawCopy("psame");
1207
1208 TH1F *hRatioS = (TH1F*)hPythia[i]->Clone(Form("hPythia_%d",i));
1209 hRatioS->Divide(hPythia[nPSmear-1]);
1210 hRatioS->SetMaximum(1E5);
1211 hRatioS->SetMinimum(0.1);
1212 hRatioS->SetYTitle("Ratio");
1213 c2->cd();
1214 if(i==0)hRatioS->DrawCopy("p");
1215 else hRatioS->DrawCopy("psame");
1216 c1->Update();
1217 c2->Update();
1218 if(getchar()=='q')return 1;
1219 }
1220
1221 if(getchar()=='q')return 1;
1222 picName = Form("jetSmearRatio");
1223 c2->SaveAs(Form(cPrintMask,picName.Data()));
1224 picName = Form("jetSmear");
1225 c1->SaveAs(Form(cPrintMask,picName.Data()));
1226 if(getchar()=='q')return 1;
1227
1228 // fetch the jet spectrum first
1229
1230
1231 c1->cd();
1232 TH1F *hPythiaRef = GetPythia8Spectrum(pName.Data(),"~/alice/sim/pythia8/examples/output/LHCJets_2.75TeV_allpt_R04.root",1.0,1);
1233 hPythiaRef->Draw();
1234 c1->Update();
1235 Float_t fMinPt[kMaxFiles];
1236 for(int iF = 0; iF < kMaxFiles;iF++){
1237 if(nColl[iF]<=0)continue;
1238 for(int i = hPythiaRef->GetNbinsX();i >0;i--){
1239 if(hPythiaRef->GetBinContent(i)>1./nColl[iF]){
1240 fMinPt[iF] = hPythiaRef->GetBinCenter(i);
1241 break;
1242 }
1243 }
1244 Printf("Ncoll %4.1f Min Pt %3.1f",nColl[iF],fMinPt[iF]);
1245 }
1246
1247 char *cFile = "~/alice/sim/pythia8/examples/output/LHCJets_2.75TeV_allpt_R04.root";
1248 TFile *fIn = (TFile*)gROOT->GetListOfFiles()->FindObject(cFile);
1249 if(!fIn) fIn = TFile::Open(cFile);
1250 if(!fIn)return 0;
1251
1252 TH1D* hPythia2[nPSmear];
1253 for(int i = 0;i < nPSmear-1;i++){
1254
1255 TH2F *hCorrelation = 0;
1256 for(int ij = 0;ij<2;ij++){
1257 TH2F* hTmp = (TH2F*)fIn->Get(Form("h2SmearCorrelationCh_C%02d_j%d",i,ij));
1258 if(hTmp&&!hCorrelation)hCorrelation = (TH2F*)hTmp->Clone(Form("%s_%d",hTmp->GetName(),i));
1259 else if (hTmp) hCorrelation->Add(hTmp);
1260 }
1261 c1->cd();
1262
1263 hPythia2[i] = hCorrelation->ProjectionY(Form("hPythia2_%d",i),10,300);
1264 hPythia2[i]->SetMarkerStyle(kFullCircle);
1265 hPythia2[i]->SetMarkerColor(kRed-9+i);
1266 hPythia2[i]->Rebin(iRebin);
1267 hPythia2[i]->SetAxisRange(0,100);
1268 hPythia2[i]->Scale(1./hPythia2[i]->GetBinWidth(1));
1269 // hist->Scale(1./xsec);
1270 hPythia2[i]->Scale(1./1.); // deta
1271
1272 // Double_t xtotal = 71.39; // xtotal for 7 TeV!!
1273 Double_t xtotal = 62.01; // xtotal for 2.76 TeV!!
1274 // scale with fraction of total NSD cross section
1275 hPythia2[i]->Scale(1/xtotal);
1276
1277
1278
1279 c1->cd();
1280 if(i==0)hPythia2[i]->DrawCopy("p");
1281 else hPythia2[i]->DrawCopy("psame");
1282
1283 TH1F *hRatioS = (TH1F*)hPythia2[i]->Clone(Form("hPythia2_%d",i));
1284 hRatioS->Divide(hPythia[nPSmear-1]);
1285 hRatioS->SetMaximum(1E5);
1286 hRatioS->SetMinimum(0.1);
1287 hRatioS->SetYTitle("Ratio");
1288 c2->cd();
1289 if(i==0)hRatioS->DrawCopy("p");
1290 else hRatioS->DrawCopy("psame");
1291 c1->Update();
1292 c2->Update();
1293 if(getchar()=='q')return 1;
1294 }
1295 return 0;
1296}
1297
1298
1299
1300void set_plot_style() {
1301 const Int_t NRGBs = 5;
1302 const Int_t NCont = 255;
1303
1304 Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
1305 Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
1306 Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
1307 Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
1308 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
1309 gStyle->SetNumberContours(NCont);
8f206bf0 1310
1311 // TGaxis::SetMaxDigits(3) // set the maximum number of axis digits (exp otherwise)
1312 // gStyle->SetPalette(1);
1313 gStyle->SetCanvasColor(10);
1314 gStyle->SetHistFillColor(10);
1315 gStyle->SetHistFillStyle(0);
1316 gStyle->SetOptTitle(0);
1317 gStyle->SetOptStat(0);
1318 gStyle->SetPadLeftMargin(0.17);
1319 gStyle->SetPadRightMargin(0.14);
1320 gStyle->SetPadBottomMargin(0.15);
1321 gStyle->SetPadTickX(1);
1322 gStyle->SetPadTickY(1);
1323 gStyle->SetAxisColor(1, "X");
1324 gStyle->SetAxisColor(1, "Y");
1325 gStyle->SetAxisColor(1, "Z");
1326 gStyle->SetLabelColor(1, "X");
1327 gStyle->SetLabelColor(1, "Y");
1328 gStyle->SetLabelColor(1, "Z");
1329 gStyle->SetTickLength(0.03, "X");
1330 gStyle->SetTickLength(0.03, "Y");
1331 gStyle->SetTickLength(0.03, "Z");
1332 gStyle->SetTitleXSize(0.07);
1333 gStyle->SetTitleYSize(0.07);
1334 gStyle->SetNdivisions(505, "X");
1335 gStyle->SetNdivisions(505, "Y");
1336 gStyle->SetNdivisions(505, "Z");
1337 gStyle->SetTitleXOffset(1.0);
1338 gStyle->SetTitleYOffset(1.0);
1339 gStyle->SetLabelOffset(0.02, "X");
1340 gStyle->SetLabelOffset(0.02, "Y");
1341 gStyle->SetLabelOffset(0.02, "Z");
1342 gStyle->SetLabelSize(0.05, "X");
1343 gStyle->SetLabelSize(0.05, "Y");
1344 gStyle->SetLabelSize(0.05, "Z");
1345 gROOT->ForceStyle();
1346
1347
e2081c22 1348}
1349
1350
1351Int_t PlotJetBFluctuations(){
1352 // plot the diffent background estimates
1353
1354 const int nCen = 4;
0282fabc 1355 const Float_t nColl[nCen] = {1502,742.5,250,46.6};
1356 const Float_t fCentLo[nCen] = {0,10,30,50};
e2081c22 1357 const Float_t fCentUp[nCen] = {10,30,50,80};
1354c5f3 1358 TH2F *hFrame = new TH2F("hFrame",";#deltap_{T} (GeV/c);Probability/GeV",200,-70,70,100,1E-5,50);
e2081c22 1359 hFrame->SetTitleOffset(1.5,"Y");
1360 hFrame->SetTitleOffset(1.5,"X");
1361 hFrame->SetLabelSize(hFrame->GetLabelSize("Y")*0.9,"Y");
1362 hFrame->SetLabelSize(hFrame->GetLabelSize("X")*0.9,"X");
1363 hFrame->SetTitleSize(hFrame->GetTitleSize("Y")*0.7,"Y");
1364 hFrame->SetTitleSize(hFrame->GetTitleSize("X")*0.7,"X");
1365
1366
1367 TString printType = "png";
1368 TString tmpName;
1369
1370 TCanvas *c1 = new TCanvas("c11","c1",600,600);
1371 c1->SetLogy();
1372
407d06f5 1373 // TFile::SetCacheFileDir("/tmp/");
e2081c22 1374
1375 // Martha, single particle jets
407d06f5 1376 TFile *fM = TFile::Open("~/alice/jets/macros/corrections/tmp/MV_PWG4_JetTasksOutput_AOD_EmbeddingSingleTrack.root");
e2081c22 1377 TH1D *hDeltaPtM[nCen] = {0};
407d06f5 1378 TString sDeltaPtM = "";
e2081c22 1379
1380 // select
407d06f5 1381 Float_t fMinPtM = 20;
1382 Float_t fMaxPtM = 40;
e2081c22 1383 int iB = 2;
1384
1385 /*
1386 0: 0-10%
1387 1: 10-30%
1388 2: 30-50%
1389 3: 50-80%
1390 */
1391
e2081c22 1392 for(int ic = 0;ic < nCen;ic++){
1393 tmpName = Form("PWG4_BkgFluctCent%dB%d",ic,iB);
1354c5f3 1394 sDeltaPtM = Form("tracks %2.0f-%2.0f GeV (MV)",fMinPtM,fMaxPtM);
e2081c22 1395 TDirectory *dir = (TDirectory*)fM->Get(tmpName.Data());
1396 if(!dir)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1397 tmpName = Form("taskBkgFluctCent%dB%d",ic,iB);
1398 TList *list = (TList*)dir->Get(tmpName.Data());
1399 if(!list)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1400
1401 TH3F *h3Tmp = (TH3F*)list->FindObject("fh3PtDeltaPtArea");
1402 hDeltaPtM[ic] = h3Tmp->ProjectionY(Form("hDeltaM%d",ic),h3Tmp->GetXaxis()->FindBin(fMinPtM),h3Tmp->GetXaxis()->FindBin(fMaxPtM),
1403 0,-1,"E");
1404 hDeltaPtM[ic]->Rebin(10);
1405 Float_t fScale = hDeltaPtM[ic]->Integral("width");
1406 if(fScale)hDeltaPtM[ic]->Scale(1./fScale);
407d06f5 1407 hDeltaPtM[ic]->SetMarkerStyle(33);
e2081c22 1408 hDeltaPtM[ic]->SetMarkerColor(kGreen+2);
1409 hDeltaPtM[ic]->SetLineColor( hDeltaPtM[ic]->GetMarkerColor());
1410 }
1411
1412 // fetch the BiAs
1413
1414 TH1D *hBiaL[nCen];
1415
407d06f5 1416 TFile *fL = TFile::Open("~/alice/jets/macros/corrections/tmp/2011-03-20_lcm_pw4plots.root");
e2081c22 1417
1418 for(int ic = 0;ic < nCen;ic++){
1419 if(iB==1)tmpName = "BiA sa antikt centrality";
1420 else if (iB==2)tmpName = "BiA va antikt centrality";
1421 TH2F *h2Tmp = (TH2F*)fL->Get(tmpName.Data());
1422 if(!h2Tmp)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1423 Int_t ibLo = h2Tmp->GetXaxis()->FindBin(fCentLo[ic]);
1424 Int_t ibUp = h2Tmp->GetXaxis()->FindBin(fCentUp[ic])-1;
1425 Printf("Line:%d bin %d - %d",__LINE__,ibLo,ibUp);
1426 hBiaL[ic] = h2Tmp->ProjectionY(Form("hBiaL%d",ic),ibLo,ibUp,"E");
1427 Float_t fScale = hBiaL[ic]->Integral("width");
1428 if(fScale) hBiaL[ic]->Scale(1./fScale);
1429 hBiaL[ic]->SetMarkerStyle(kFullCircle);
1430 hBiaL[ic]->SetMarkerColor(kBlue);
1431 hBiaL[ic]->SetLineColor( hBiaL[ic]->GetMarkerColor());
1432 }
1433
1434 TH1D *hBiaRC[nCen];
407d06f5 1435 TString sBiaRC = "";
e2081c22 1436
1437 for(int ic = 0;ic < nCen;ic++){
407d06f5 1438 if(iB==1)tmpName = "BiA sa RC skip0 centrality";
1439 else if (iB==2)tmpName = "BiA va RC skip0 centrality";
1440 sBiaRC = "BiA RC";
e2081c22 1441 TH2F *h2Tmp = (TH2F*)fL->Get(tmpName.Data());
1442 if(!h2Tmp)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1443 Int_t ibLo = h2Tmp->GetXaxis()->FindBin(fCentLo[ic]);
1444 Int_t ibUp = h2Tmp->GetXaxis()->FindBin(fCentUp[ic])-1;
1445 Printf("Line:%d bin %d - %d",__LINE__,ibLo,ibUp);
1446 hBiaRC[ic] = h2Tmp->ProjectionY(Form("hBiaRC%d",ic),ibLo,ibUp,"E");
1447 Float_t fScale = hBiaRC[ic]->Integral("width");
1448 if(fScale) hBiaRC[ic]->Scale(1./fScale);
1449 hBiaRC[ic]->SetMarkerStyle(kFullCircle);
1450 hBiaRC[ic]->SetMarkerColor(kRed);
1451 hBiaRC[ic]->SetLineColor( hBiaRC[ic]->GetMarkerColor());
1452 }
1453
1454
407d06f5 1455 TH1D *hBiaRC2[nCen];
1456 TString sBiaRC2;
1457 for(int ic = 0;ic < nCen;ic++){
1458 if(iB==1)tmpName = "BiA sa RC skip2 centrality";
1459 else if (iB==2)tmpName = "BiA va RC skip2 centrality";
1460 sBiaRC2 = "BiA RC (excl. 2 leading jets)";
1461 TH2F *h2Tmp = (TH2F*)fL->Get(tmpName.Data());
1462 if(!h2Tmp)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1463 Int_t ibLo = h2Tmp->GetXaxis()->FindBin(fCentLo[ic]);
1464 Int_t ibUp = h2Tmp->GetXaxis()->FindBin(fCentUp[ic])-1;
1465 Printf("Line:%d bin %d - %d",__LINE__,ibLo,ibUp);
1466 hBiaRC2[ic] = h2Tmp->ProjectionY(Form("hBiaRC2%d",ic),ibLo,ibUp,"E");
1467 Float_t fScale = hBiaRC2[ic]->Integral("width");
1468 if(fScale) hBiaRC2[ic]->Scale(1./fScale);
1469 hBiaRC2[ic]->SetMarkerStyle(kOpenCircle);
1470 hBiaRC2[ic]->SetMarkerColor(kRed);
1471 hBiaRC2[ic]->SetLineColor( hBiaRC2[ic]->GetMarkerColor());
1472 }
1473
1474
e2081c22 1475
1476 c1->SetLogy(0);
1477 c1->SetLogz();
1478 TLatex *txt = new TLatex();
1479 txt->SetTextFont(gStyle->GetTextFont());
1480 txt->SetTextSize(gStyle->GetTextSize()*0.6);
1481
1482 TLatex *txt2 = new TLatex();
1483 txt2->SetTextFont(gStyle->GetTextFont());
1484 txt2->SetTextSize(gStyle->GetTextSize()*0.7);
1485 txt2->SetTextAlign(22);
1486 txt2->SetTextColor(kRed);
1487
1488
1489 if(iB==1)tmpName = "background sa vs multiplicity";
1490 else if(iB==2)tmpName = "background va vs multiplicity";
1491 TH2F *h2RhoVsMult = (TH2F*)fL->Get(tmpName.Data());
1492 if(!h2RhoVsMult)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1493 c1->SetMargin(0.15,0.18,0.2,0.10);
1494
1495 h2RhoVsMult->SetTitleOffset(1.5,"Y");
1496 h2RhoVsMult->SetTitleOffset(1.5,"X");
1497 h2RhoVsMult->SetLabelSize(h2RhoVsMult->GetLabelSize("Y")*0.9,"Y");
1498 h2RhoVsMult->SetLabelSize(h2RhoVsMult->GetLabelSize("X")*0.9,"X");
1499 h2RhoVsMult->SetTitleSize(h2RhoVsMult->GetTitleSize("Y")*0.7,"Y");
1500 h2RhoVsMult->SetTitleSize(h2RhoVsMult->GetTitleSize("X")*0.7,"X");
1501 h2RhoVsMult->SetXTitle("input tracks");
1502 h2RhoVsMult->SetYTitle("#rho (GeV/unit area)");
1503 h2RhoVsMult->SetAxisRange(0,3000.);
1504 h2RhoVsMult->Draw("colz");
1505 txt->Draw();
1354c5f3 1506 txt->DrawLatex(100,180,"LHC 2010 Pb-Pb Run #sqrt{s_{NN}} = 2.76 TeV");
407d06f5 1507 // txt2->DrawLatex(800,150,"ALICE Performance");
1508 // txt2->DrawLatex(800,140,"01/03/2011");
e2081c22 1509 c1->Update();
1510 c1->SaveAs(Form("rhovsmult_B%d.%s",iB,printType.Data()));
1511 if(getchar()=='q')return 1;
1512
1513 if(iB==1)tmpName = "background sa vs centrality";
1514 else if(iB==2)tmpName = "background va vs centrality";
1515 TH2F *h2RhoVsCent = (TH2F*)fL->Get(tmpName.Data());
1516 if(!h2RhoVsCent)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1517 h2RhoVsCent->SetXTitle("centrality (%)");
1518 h2RhoVsCent->SetYTitle("#rho (GeV/unit area)");
1519
1520 h2RhoVsCent->SetTitleOffset(1.5,"Y");
1521 h2RhoVsCent->SetTitleOffset(1.5,"X");
1522 h2RhoVsCent->SetLabelSize(h2RhoVsCent->GetLabelSize("Y")*0.9,"Y");
1523 h2RhoVsCent->SetLabelSize(h2RhoVsCent->GetLabelSize("X")*0.9,"X");
1524 h2RhoVsCent->SetTitleSize(h2RhoVsCent->GetTitleSize("Y")*0.7,"Y");
1525 h2RhoVsCent->SetTitleSize(h2RhoVsCent->GetTitleSize("X")*0.7,"X");
1526 h2RhoVsCent->SetAxisRange(3,200,"Y");
1527 h2RhoVsCent->Draw("colz");
1354c5f3 1528 // txt->DrawLatex(20,180,"LHC 2010 Pb-Pb Run #sqrt{s_{NN}} = 2.76 TeV");
407d06f5 1529 // txt2->DrawLatex(50,150,"ALICE Performance");
1530 // txt2->DrawLatex(50,140,"01/03/2011");
e2081c22 1531 c1->Update();
1532 c1->SaveAs(Form("rhovscent_B%d.%s",iB,printType.Data()));
1533 if(getchar()=='q')return 1;
1534
1535 // fetch the data from bastian...
1536 Float_t fMinPtB = fMinPtM;
1537 Float_t fMaxPtB = fMaxPtM;
1538
1539
e2081c22 1540 TH1D *hDeltaPtB1[nCen] = {0};
407d06f5 1541 TString sDeltaPtB1 = "";
e2081c22 1542 for(int ic = 0;ic < nCen;ic++){
407d06f5 1543 sDeltaPtB1 = Form("anti-k_{T} embedded jet %2.0f-%2.0f GeV",fMinPtB,fMaxPtB);
8dc333da 1544 TH2F *hTmp = GetTH2PlotB("~/Dropbox/SharedJets/Bastian/Files/110428/",1,0,ic); // emb jets
e2081c22 1545 if(!hTmp)Printf("%d %s not found",__LINE__,tmpName.Data());
1546 int ibLo = hTmp->GetYaxis()->FindBin(fMinPtB);
1547 int ibUp = hTmp->GetYaxis()->FindBin(fMaxPtB)-1;
1548 hDeltaPtB1[ic] = hTmp->ProjectionX(Form("fHistDeltaPtB1_c%d",ic),ibLo,ibUp,"E");
1549 hDeltaPtB1[ic]->SetMarkerStyle(kFullSquare);
1550 hDeltaPtB1[ic]->SetMarkerColor(kBlue);
1551 hDeltaPtB1[ic]->SetLineColor(hDeltaPtB1[ic]->GetMarkerColor());
1552 hDeltaPtB1[ic]->Rebin(2);
1553 Float_t fScale = hDeltaPtB1[ic]->Integral("width");
1554 if(fScale) hDeltaPtB1[ic]->Scale(1./fScale);
1555 }
1556
1557
407d06f5 1558 TH1D *hDeltaPtB2[nCen] = {0};
1559 TString sDeltaPtB2 = "";
1560 for(int ic = 0;ic < nCen;ic++){
8f206bf0 1561 sDeltaPtB2 = Form("anti-k_{T} embedded track %2.0f-%2.0f GeV",fMinPtB,fMaxPtB);
8dc333da 1562 TH2F *hTmp = GetTH2PlotB("~/Dropbox/SharedJets/Bastian/Files/110430/",0,0,ic); // emb jets
407d06f5 1563 if(!hTmp)Printf("%d %s not found",__LINE__,tmpName.Data());
1564 int ibLo = hTmp->GetYaxis()->FindBin(fMinPtB);
1565 int ibUp = hTmp->GetYaxis()->FindBin(fMaxPtB)-1;
1566 hDeltaPtB2[ic] = hTmp->ProjectionX(Form("fHistDeltaPtB2_c%d",ic),ibLo,ibUp,"E");
1567 hDeltaPtB2[ic]->SetMarkerStyle(33);
1568 hDeltaPtB2[ic]->SetMarkerColor(kBlue+4);
1569 hDeltaPtB2[ic]->SetLineColor(hDeltaPtB2[ic]->GetMarkerColor());
1570 hDeltaPtB2[ic]->Rebin(2);
1571 Float_t fScale = hDeltaPtB2[ic]->Integral("width");
1572 if(fScale) hDeltaPtB2[ic]->Scale(1./fScale);
1573 }
1574
1575
1576
e2081c22 1577 c1->SetLogy();
1578 c1->SetMargin(0.15,0.05,0.2,0.05);
1579
1580 TF1 *gaus = new TF1("gaus","gaus",-60,2);
1581 TF1 *gaus2 = new TF1("gaus2","gaus",-60,2);
407d06f5 1582 Double_t mean = 0;
1583 Double_t sigma = 0;
1584 Double_t sigma_err = 0;
1585 TF1* tmpGaus = 0;
c08c3ad2 1586 TF1 *gamma0 = new TF1("gamma0",Gamma0,-40,40,4);
1587 gamma0->SetParameters(1,1,1,1);
407d06f5 1588
1354c5f3 1589
1590
e2081c22 1591 for(int ic = 0;ic < nCen;ic++){
407d06f5 1592 TLegend *leg1 = new TLegend(0.2,0.65,0.3,0.93);
1354c5f3 1593
1594 leg1->SetHeader(Form("Pb-Pb %2.0f-%2.0f%% R = 0.4 (B%d)",fCentLo[ic],fCentUp[ic],iB));
e2081c22 1595 leg1->SetFillColor(0);
1596 leg1->SetTextFont(gStyle->GetTextFont());
1597 leg1->SetTextSize(gStyle->GetTextSize()*0.6);
1598 leg1->SetBorderSize(0);
1599 hFrame->DrawCopy();
1600
1601 /*
1602 hBiaL[ic]->DrawCopy("psame");
1603 leg1->AddEntry(hBiaL[ic],Form("BiA anti-k_{T}"),"P");
1604 */
1605
c08c3ad2 1606 /*
407d06f5 1607 hBiaRC2[ic]->DrawCopy("psame");
1608 tmpGaus = FitLHSgaus(hBiaRC2[ic],sigma,sigma_err,mean);
1609 tmpGaus->SetRange(-40,40);
1610 tmpGaus->SetLineColor( hBiaRC2[ic]->GetLineColor());
1611 tmpGaus->SetLineStyle(kDashed);
1612 tmpGaus->Draw("same");
1613 leg1->AddEntry(hBiaRC2[ic],sBiaRC2.Data(),"P");
1614 leg1->AddEntry(tmpGaus,Form("LHS Gaus fit: #mu = %2.2f, #sigma = %2.2f",mean,sigma),"L");
c08c3ad2 1615 */
1616
e2081c22 1617 hBiaRC[ic]->DrawCopy("psame");
c08c3ad2 1618 hBiaRC[ic]->Fit(gamma0);
1619 c1->Update();
1620 if(getchar()=='q')return 1;
1621 continue;
1622
1623
407d06f5 1624 tmpGaus = FitLHSgaus(hBiaRC[ic],sigma,sigma_err,mean);
1625 tmpGaus->SetRange(-40,40);
1626 tmpGaus->SetLineColor( hBiaRC[ic]->GetLineColor());
1627 tmpGaus->Draw("same");
1628 leg1->AddEntry(hBiaRC[ic],sBiaRC.Data(),"P");
1354c5f3 1629 leg1->AddEntry(tmpGaus,Form("LHS Gauss fit: #mu = %2.2f, #sigma = %2.2f",mean,sigma),"L");
407d06f5 1630
1631 hDeltaPtB1[ic]->Draw("psame");
1632 tmpGaus = FitLHSgaus(hDeltaPtB1[ic],sigma,sigma_err,mean);
1633 tmpGaus->SetRange(-40,40);
1634 tmpGaus->SetLineColor( hDeltaPtB1[ic]->GetLineColor());
1635 tmpGaus->Draw("same");
1636 leg1->AddEntry(hDeltaPtB1[ic],sDeltaPtB1.Data(),"P");
1354c5f3 1637 leg1->AddEntry(tmpGaus,Form("LHS Gauss fit: #mu = %2.2f, #sigma = %2.2f",mean,sigma),"L");
407d06f5 1638
1639
1640 hDeltaPtB2[ic]->Draw("psame");
1641 tmpGaus = FitLHSgaus(hDeltaPtB2[ic],sigma,sigma_err,mean);
1642 tmpGaus->SetRange(-40,40);
1643 tmpGaus->SetLineColor( hDeltaPtB2[ic]->GetLineColor());
1644 tmpGaus->Draw("same");
1645 leg1->AddEntry(hDeltaPtB2[ic],sDeltaPtB2.Data(),"P");
1354c5f3 1646 leg1->AddEntry(tmpGaus,Form("LHS Gauss fit: #mu = %2.2f, #sigma = %2.2f",mean,sigma),"L");
e2081c22 1647
1648 hDeltaPtM[ic]->DrawCopy("psame");
407d06f5 1649 tmpGaus = FitLHSgaus(hDeltaPtM[ic],sigma,sigma_err,mean);
1650 tmpGaus->SetRange(-40,40);
1651 tmpGaus->SetLineColor( hDeltaPtM[ic]->GetLineColor());
1652 tmpGaus->Draw("same");
1653 leg1->AddEntry(hDeltaPtM[ic],sDeltaPtM.Data(),"P");
1354c5f3 1654 leg1->AddEntry(tmpGaus,Form("LHS Gauss fit: #mu = %2.2f, #sigma = %2.2f",mean,sigma),"L");
407d06f5 1655
e2081c22 1656 leg1->Draw();
407d06f5 1657
1658
e2081c22 1659 txt2->SetNDC();
407d06f5 1660 // txt2->DrawLatex(0.5,0.97,"ALICE Performance 01/03/2011");
e2081c22 1661 c1->Update();
407d06f5 1662 c1->SaveAs(Form("deltaPt_pT_%d_%d_B%d_cen%02d.%s",(Int_t)fMinPtM,(Int_t)fMaxPtM,iB,ic,printType.Data()));
e2081c22 1663 if(getchar()=='q')return 1;
1664
1665
1666 }
1667
1668
1669
0282fabc 1670 // now plot different centralities on one plot
1671 bool bFirst = true;
1672 hFrame->DrawCopy();
1673 TLegend *leg2 = new TLegend(0.2,0.65,0.3,0.93);
1674 leg2->SetHeader(Form("%s R = 0.4 (B%d)",sBiaRC2.Data(),iB));
1675 TH1D **hTmp1 = hBiaRC2;
8dc333da 1676 tmpName = "BiARC2";
0282fabc 1677 leg2->SetFillColor(0);
1678 leg2->SetTextFont(gStyle->GetTextFont());
1679 leg2->SetTextSize(gStyle->GetTextSize()*0.6);
1680 leg2->SetBorderSize(0);
1681 hFrame->DrawCopy();
1682
1683
1684 for(int ic = 0;ic<nCen;ic++){
1685 int iColor = hTmp1[ic]->GetMarkerColor()-4+2*ic;
1686 hTmp1[ic]->SetLineColor(iColor);
1687 hTmp1[ic]->SetMarkerColor(iColor);
1688 if(bFirst){
1689 hTmp1[ic]->Draw("psame");
1690 bFirst = false;
8dc333da 1691 }
0282fabc 1692 else{
1693 hTmp1[ic]->Draw("psame");
1694 }
1695 leg2->AddEntry(hTmp1[ic],Form("%2.0f-%2.0f%%",fCentLo[ic],fCentUp[ic]),"P");
1696
1697 }
1698 leg2->Draw("same");
1699 c1->Update();
1700 c1->SaveAs(Form("deltaPt_cent_%s_B%d.%s",tmpName.Data(),iB,printType.Data()));
8dc333da 1701 if(getchar()=='q')return 1;
0282fabc 1702
1703 bFirst = true;
1704 for(int ic = nCen-1;ic>=0;ic--){
1705 hTmp1[ic]->Scale(1./nColl[ic]);
1706 hTmp1[ic]->SetMinimum(1E-8);
1707 if(bFirst){
1708 hTmp1[ic]->Draw("p");
1709 bFirst = false;
8dc333da 1710 }
0282fabc 1711 else{
1712 hTmp1[ic]->Draw("psame");
1713 }
1714
1715 }
1716 leg2->Draw("same");
1717 c1->Update();
1718 c1->SaveAs(Form("deltaPt_cent_ncoll_%s_B%d.%s",tmpName.Data(),iB,printType.Data()));
1719
1720
1721
e2081c22 1722 return 0;
1723
1724
1725}
1726
8dc333da 1727
1354c5f3 1728Int_t PlotJetBFluctuations2(UInt_t iPlotFlag,UInt_t iPlotType,Float_t inPtLo,Float_t inPtUp,int iALICEType){
8dc333da 1729 // plot the diffent background estimates
1730
1731 const int nCen = 4;
1732 const Float_t nColl[nCen] = {1502,742.5,250,46.6};
1733 const Float_t fCentLo[nCen] = {0,10,30,50};
1734 const Float_t fCentUp[nCen] = {10,30,50,80};
1735
1736 Int_t iB = 2;
1737
1738
1739
1740 // multiplicity (nb. of tracklets) bins
1741 const Int_t nMult = 17;
1742 const Int_t multBinWidth = 200;
1743 Int_t multMin[nMult] = {0,};
1744 Int_t multMax[nMult] = {0,};
1745 for(Int_t i=0; i<nMult; ++i){
1746 multMin[i] = i*multBinWidth;
1747 multMax[i] = multMin[i]+multBinWidth-1;
1748 }
1749
1750 // bins wrt. reactions
1751 const Int_t nRP = 4; // 0 --> all 1,2,3
1752
1354c5f3 1753 TH2F *hFrame = new TH2F("hFrame",";#deltap_{T} (GeV/c);probability density (c/GeV)",200,-70,100,1000,1E-7,50);
8dc333da 1754 hFrame->SetTitleOffset(1.5,"Y");
1755 hFrame->SetTitleOffset(1.5,"X");
1756 hFrame->SetLabelSize(hFrame->GetLabelSize("Y")*0.9,"Y");
1757 hFrame->SetLabelSize(hFrame->GetLabelSize("X")*0.9,"X");
1758 hFrame->SetTitleSize(hFrame->GetTitleSize("Y")*0.7,"Y");
1759 hFrame->SetTitleSize(hFrame->GetTitleSize("X")*0.7,"X");
1760
8f206bf0 1761 TCanvas *c1 = new TCanvas("c11","c1",800,800);
8dc333da 1762 c1->SetLogy();
1763
1764
1765 // Produce Delta Pt Plots for each centrality and for each type
5f5e3e77 1766 const int kDeltaTypes = 12;
8dc333da 1767 TString sDelta[kDeltaTypes][nRP];
1768 TH1D *hDeltaPt[kDeltaTypes][nCen][nRP] = {0,};
1769 TGraphErrors *grMeanDeltaPtCent[kDeltaTypes][nRP] = {0,};
1770 TGraphErrors *grSigmaDeltaPtCent[kDeltaTypes][nRP] = {0,};
1771
1772 //
1773 TH1D *hDeltaPtMult[kDeltaTypes][nMult][nRP] = {0,};
1774 TGraphErrors *grMeanDeltaPtMult[kDeltaTypes][nRP] = {0,};
1775 TGraphErrors *grSigmaDeltaPtMult[kDeltaTypes][nRP] = {0,};
1776
1777 //
5f5e3e77 1778 const Int_t kMarkerDelta[kDeltaTypes] = {kFullCircle,kFullCircle,kOpenCircle, 33,kFullSquare,kFullSquare,kFullSquare, kFullSquare, 31,kFullCircle,kOpenCircle,31};
1354c5f3 1779 const Int_t kColorDelta[kDeltaTypes] = {kRed, kBlue, kRed, kBlue, kBlue+2, kGray+1, kRed+2, kOrange+2,kOrange+2,kRed,kRed,kBlue};
5f5e3e77 1780 const Int_t kFitDelta[kDeltaTypes] = { 1, 1, 1, 1, 1, 1, 1, 0, 1,1,1,0};
1781 const Int_t kNormDelta[kDeltaTypes] = { 2, 0, 0, 2, 2, 0, 0, 0, 1,2,0,0};
1354c5f3 1782 const Int_t kLineDelta[kDeltaTypes] = { 1, 1, 2, 1, 1, 1, 1, 1, 1,1,1,1};
8dc333da 1783 const Int_t kMarkerRP[nRP] = {kFullSquare,22,29,23}; // first no used
1784 const Int_t kColorRP[nRP] = {kBlack,kBlue+2,kCyan+2,kGreen+2};
1785
1786 //
1787
1788
1789 const Int_t kNormType = 7; // normaize this spectrum at fixed p_T
8f206bf0 1790 const Float_t kNormPt = 20.; // normaize this spectrum at fixed p_T
8dc333da 1791 Float_t kNormValue[nMult][nCen] = {0,}; // normaize the spectrum with this value, if negative take from the first delta p_T histo
1792 Float_t kNormValueMult[nMult][nRP] = {0,}; // normaize the spectrum with this value, if negative take from the first delta p_T histo
1793
1794
1795 //
1796 // the cuts for the single centrality
1797 Float_t fMinPt[kDeltaTypes] = {0};
1798 Float_t fMaxPt[kDeltaTypes] = {0};
1799
1800 for(int iD = 0;iD<kDeltaTypes;iD++){
1801 fMinPt[iD] = inPtLo;
1802 fMaxPt[iD] = inPtUp;
1803 }
1804
1805 TString tmpName;
1806
1807
1808
1809 // 0 Leticias BiA anti-kT
1810
1811 TFile *fL = TFile::Open("~/Dropbox/SharedJets/Leticia/randomcones/pwg4plots.root");
1812 TFile *fRP = TFile::Open("~/Dropbox/SharedJets/Leticia/randomcones/reactionplane.root");
5f5e3e77 1813 TFile *fMRP = TFile::Open("~/Dropbox/SharedJets/Leticia/randomconesvsmult.root");
8dc333da 1814 if(fDebug)Printf("Line: %d",__LINE__);
1815 Int_t iDelta = 0;
1816 Int_t iRP = 0;
5f5e3e77 1817 Float_t fScale = 1;
1818
8dc333da 1819 if(fL){
1820 iDelta = 0;
1821 // leticia BiA Randome cones
1354c5f3 1822 sDelta[iDelta][0] = "rand. cones";
8dc333da 1823 if(iB==1)tmpName = "BiA sa RC skip0 centrality";
1824 else if (iB==2)tmpName = "BiA RC skip0 va centrality";
1825 TH2F *h2Tmp = (TH2F*)fL->Get(tmpName.Data());
1826 if(!h2Tmp)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1827
1828
1829 for(int ic = 0;ic < nCen;ic++){
1830 iRP = 0;
1831 if(!h2Tmp)continue;
1832 Int_t ibLo = h2Tmp->GetXaxis()->FindBin(fCentLo[ic]);
1833 Int_t ibUp = h2Tmp->GetXaxis()->FindBin(fCentUp[ic])-1;
1834 hDeltaPt[iDelta][ic][iRP] = h2Tmp->ProjectionY(Form("hBiaRC%d",ic),ibLo,ibUp,"E");
5f5e3e77 1835
1836
1837 fScale = hDeltaPt[iDelta][ic][iRP]->Integral("width");
8dc333da 1838 if(fScale) hDeltaPt[iDelta][ic][iRP]->Scale(1./fScale);
1839 if(fRP){
1840 if(ic==0){tmpName = "bia RC vs RP central";}
1841 else if(ic == 3){ tmpName = "bia RC vs RP perif";}
1842 else {continue;}
1843 TH2F *h2TmpRP = (TH2F*)fRP->Get(tmpName.Data());
1844 if(!h2TmpRP)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1845 for(iRP = 1;iRP<nRP;iRP++){
1846 hDeltaPt[iDelta][ic][iRP] = h2TmpRP->ProjectionY(Form("hDeltaPt%d_%d_%d",iDelta,ic,iRP),iRP,iRP);
5f5e3e77 1847 if(kNormDelta[iDelta]==2){
1848
1849 }
1850 else{
1851 fScale = hDeltaPt[iDelta][ic][iRP]->Integral("width");
1852 }
8dc333da 1853 if(fScale) hDeltaPt[iDelta][ic][iRP]->Scale(1./fScale);
1854 }
1855 }
1856 }
5f5e3e77 1857
8dc333da 1858 iRP = 0;
1859
1860
1861 if(fDebug)Printf("Line: %d",__LINE__);
5f5e3e77 1862
8dc333da 1863
8dc333da 1864 for(int im = 0;im < nMult;im++){
5f5e3e77 1865 for(iRP=1;iRP<nRP;iRP++){
1866 tmpName = Form("hBiaRP%d_M%d_cut0.15",iRP-1,im);
1867 hDeltaPtMult[iDelta][im][iRP] = (TH1D*)fMRP->Get(tmpName.Data());
1868 if(!hDeltaPtMult[iDelta][im][iRP] )continue;
1869 if(iRP==1)hDeltaPtMult[iDelta][im][0] = (TH1D*)hDeltaPtMult[iDelta][im][iRP]->Clone( Form("hBiaRP%d_M%d_cut0.15",-1,im));
1870 else hDeltaPtMult[iDelta][im][0]->Add(hDeltaPtMult[iDelta][im][iRP]);
1871 }
1872 // norm
1873 for(iRP=0;iRP<nRP;iRP++){
1874 if(kNormDelta[iDelta]==2){
1875 if(iRP==0)fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
1876 }
1877 else{
1878 fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
1879 }
1880 if(fScale) hDeltaPtMult[iDelta][im][iRP]->Scale(1./fScale);
1881 }
8dc333da 1882 }
5f5e3e77 1883
1884 iRP = 0;
8dc333da 1885 if(fDebug)Printf("Line: %d",__LINE__);
1886 //------------
1887
1888 iDelta = 1;
1354c5f3 1889 sDelta[iDelta][0] = "RC (randomized #eta#phi)";
8dc333da 1890 if(iB==1)tmpName = "BiA RC random input sa centrality";
1891 else if (iB==2)tmpName = "BiA RC random input va centrality";
1892 h2Tmp = (TH2F*)fL->Get(tmpName.Data());
1893 if(!h2Tmp)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1894 if(fDebug)Printf("Line: %d",__LINE__);
1895 for(int ic = 0;ic < nCen;ic++){
1354c5f3 1896 if(!h2Tmp)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
8dc333da 1897 if(!h2Tmp)continue;
1898 iRP = 0;
1899 Int_t ibLo = h2Tmp->GetXaxis()->FindBin(fCentLo[ic]);
1900 Int_t ibUp = h2Tmp->GetXaxis()->FindBin(fCentUp[ic])-1;
1901 Printf("Line:%d bin %d - %d",__LINE__,ibLo,ibUp);
1354c5f3 1902 hDeltaPt[iDelta][ic][iRP] = h2Tmp->ProjectionY(Form("hBia1L%d",ic),ibLo,ibUp,"E");
8dc333da 1903 Float_t fScale = hDeltaPt[iDelta][ic][iRP]->Integral("width");
1904 if(fScale) hDeltaPt[iDelta][ic][iRP]->Scale(1./fScale);
1905 if(fRP){
1906 if(ic==0){tmpName = "bia RC randomized input vs RP central";}
1907 else if(ic == 3){ tmpName = "bia RC randomized input vs RP perif";}
1908 else {continue;}
1909 TH2F *h2TmpRP = (TH2F*)fRP->Get(tmpName.Data());
1910 if(!h2TmpRP)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1911 for(iRP = 1;iRP<nRP;iRP++){
1912 hDeltaPt[iDelta][ic][iRP] = h2TmpRP->ProjectionY(Form("hDeltaPt%d_%d_%d",iDelta,ic,iRP),iRP,iRP);
1913 Float_t fScale = hDeltaPt[iDelta][ic][iRP]->Integral("width");
1914 if(fScale) hDeltaPt[iDelta][ic][iRP]->Scale(1./fScale);
1915 }
1916 }
1354c5f3 1917 Printf(hDeltaPt[iDelta][ic][iRP])
8dc333da 1918 }
8dc333da 1919
1354c5f3 1920 iRP = 0;
5f5e3e77 1921 for(int im = 0;im < nMult;im++){
1922 for(iRP=1;iRP<nRP;iRP++){
1923 tmpName = Form("hBia_random_RP%d_M%d_cut0.15",iRP-1,im);
1924 hDeltaPtMult[iDelta][im][iRP] = (TH1D*)fMRP->Get(tmpName.Data());
1925 if(!hDeltaPtMult[iDelta][im][iRP] )continue;
1926 if(iRP==1)hDeltaPtMult[iDelta][im][0] = (TH1D*)hDeltaPtMult[iDelta][im][iRP]->Clone( Form("hBiaRP%d_M%d_cut0.15",-1,im));
1927 else hDeltaPtMult[iDelta][im][0]->Add(hDeltaPtMult[iDelta][im][iRP]);
1928 }
1929 // norm
1930 for(iRP=0;iRP<nRP;iRP++){
1931 if(kNormDelta[iDelta]==2){
1932 if(iRP==0)fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
1933 }
1934 else{
1935 fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
1936 }
1937 if(fScale) hDeltaPtMult[iDelta][im][iRP]->Scale(1./fScale);
1938 }
1939 }
1940
1941 /*
8dc333da 1942 if(fDebug)Printf("Line: %d",__LINE__);
1943 if(iB==1)tmpName = "BiA RC random input sa multiplicity";
1944 else if (iB==2)tmpName = "BiA RC random input va multiplicity";
1945 h2Tmp = (TH2F*)fL->Get(tmpName.Data());
1946 if(!h2Tmp)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1947 for(int im = 0;im < nMult;im++){
1948 if(!h2Tmp)continue;
1949 Int_t ibLo = h2Tmp->GetXaxis()->FindBin(multMin[im]);
1950 Int_t ibUp = h2Tmp->GetXaxis()->FindBin(multMax[im])-1;
1951 Printf("Line:%d bin %d - %d",__LINE__,ibLo,ibUp);
1952 hDeltaPtMult[iDelta][im][iRP] = h2Tmp->ProjectionY(Form("hBia%d_M%d",iDelta,im),ibLo,ibUp,"E");
1953 Float_t fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
1954 if(fScale) hDeltaPtMult[iDelta][im][iRP]->Scale(1./fScale);
1955 }
5f5e3e77 1956 */
8dc333da 1957
1958
1959 // ------------------------------------------
1960 // leticia random cones skip leading jets
1354c5f3 1961 iDelta = 2;iRP = 0;
1962 sDelta[iDelta][0] = "RC (w/o 2 lead. jets)";
8dc333da 1963 if(fDebug)Printf("Line: %d",__LINE__);
1964 if(iB==1)tmpName = "BiA sa RC skip2 centrality";
1965 else if (iB==2)tmpName = "BiA RC skip2 va centrality";
1966 h2Tmp = (TH2F*)fL->Get(tmpName.Data());
1967 if(!h2Tmp)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1968 for(int ic = 0;ic < nCen;ic++){
1969 if(!h2Tmp)continue;
8dc333da 1970 Int_t ibLo = h2Tmp->GetXaxis()->FindBin(fCentLo[ic]);
1971 Int_t ibUp = h2Tmp->GetXaxis()->FindBin(fCentUp[ic])-1;
1972 Printf("Line:%d bin %d - %d",__LINE__,ibLo,ibUp);
1973 hDeltaPt[iDelta][ic][iRP] = h2Tmp->ProjectionY(Form("hBiaRC2%d",ic),ibLo,ibUp,"E");
1974 Float_t fScale = hDeltaPt[iDelta][ic][iRP]->Integral("width");
1975 if(fScale) hDeltaPt[iDelta][ic][iRP]->Scale(1./fScale);
1976 }
5f5e3e77 1977
8dc333da 1978 if(fDebug)Printf("Line: %d",__LINE__);
1979 if(iB==1)tmpName = "BiA RC skip2 sa multiplicity";
1980 else if (iB==2)tmpName = "BiA RC skip2 va multiplicity";
1981 h2Tmp = (TH2F*)fL->Get(tmpName.Data());
1982 if(!h2Tmp)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
1983 for(int im = 0;im < nMult;im++){
1984 if(!h2Tmp)continue;
1985 Int_t ibLo = h2Tmp->GetXaxis()->FindBin(multMin[im]);
1986 Int_t ibUp = h2Tmp->GetXaxis()->FindBin(multMax[im])-1;
1987 Printf("Line:%d bin %d - %d",__LINE__,ibLo,ibUp);
1988 hDeltaPtMult[iDelta][im][iRP] = h2Tmp->ProjectionY(Form("hBia%d_M%d",iDelta,im),ibLo,ibUp,"E");
1989 Float_t fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
1990 if(fScale) hDeltaPtMult[iDelta][im][iRP]->Scale(1./fScale);
5f5e3e77 1991 }
8dc333da 1992
8dc333da 1993
5f5e3e77 1994 //
1995 iDelta = 9;
1996 sDelta[iDelta][0] = "RC (cut 2 GeV)";
1997 for(int im = 0;im < nMult;im++){
1998 for(iRP=1;iRP<nRP;iRP++){
1999 tmpName = Form("hBiaRP%d_M%d_cut2",iRP-1,im);
2000 hDeltaPtMult[iDelta][im][iRP] = (TH1D*)fMRP->Get(tmpName.Data());
2001 if(!hDeltaPtMult[iDelta][im][iRP] )continue;
2002 if(iRP==1)hDeltaPtMult[iDelta][im][0] = (TH1D*)hDeltaPtMult[iDelta][im][iRP]->Clone( Form("hBiaRP%d_M%d_cut0.15",-1,im));
2003 else hDeltaPtMult[iDelta][im][0]->Add(hDeltaPtMult[iDelta][im][iRP]);
2004 }
2005 // norm
2006 for(iRP=0;iRP<nRP;iRP++){
2007 if(kNormDelta[iDelta]==2){
2008 if(iRP==0)fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
2009 }
2010 else{
2011 fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
2012 }
2013 if(fScale) hDeltaPtMult[iDelta][im][iRP]->Scale(1./fScale);
2014 }
2015 }
8dc333da 2016
5f5e3e77 2017 iDelta = 10;
2018 sDelta[iDelta][0] = "RC randomized event (cut 2 GeV)";
2019 for(int im = 0;im < nMult;im++){
2020 for(iRP=1;iRP<nRP;iRP++){
2021 tmpName = Form("hBia_random_RP%d_M%d_cut2",iRP-1,im);
2022 hDeltaPtMult[iDelta][im][iRP] = (TH1D*)fMRP->Get(tmpName.Data());
2023 if(!hDeltaPtMult[iDelta][im][iRP] )continue;
2024 if(iRP==1)hDeltaPtMult[iDelta][im][0] = (TH1D*)hDeltaPtMult[iDelta][im][iRP]->Clone( Form("hBiaRP%d_M%d_cut0.15",-1,im));
2025 else hDeltaPtMult[iDelta][im][0]->Add(hDeltaPtMult[iDelta][im][iRP]);
2026 }
2027 // norm
2028 for(iRP=0;iRP<nRP;iRP++){
2029 if(kNormDelta[iDelta]==2){
2030 if(iRP==0)fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
2031 }
2032 else{
2033 fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
2034 }
2035 if(fScale) hDeltaPtMult[iDelta][im][iRP]->Scale(1./fScale);
2036 }
8dc333da 2037 }
5f5e3e77 2038
2039
8dc333da 2040 }
5f5e3e77 2041 // fetch the data from bastian...
2042
2043
2044 // jet embedded
2045 TString cBB = "/Users/kleinb/Dropbox/SharedJets/Bastian/Files/";
2046 iDelta = 3;
1354c5f3 2047 sDelta[iDelta][0] = Form("track: %2.0f-%2.0f GeV",fMinPt[iDelta],fMaxPt[iDelta]);
5f5e3e77 2048 normRP = 1;
2049 for(iRP = 0;iRP<nRP;iRP++){
2050 for(int ic = 0;ic < nCen;ic++){
2051 TH2F *hTmp = GetTH2PlotB(cBB.Data(),0,0,ic,iRP-1); // emb jets
2052 if(!hTmp)Printf("%d %s not found",__LINE__,tmpName.Data());
2053 int ibLo = hTmp->GetYaxis()->FindBin(fMinPt[iDelta]);
2054 int ibUp = hTmp->GetYaxis()->FindBin(fMaxPt[iDelta])-1;
2055 hDeltaPt[iDelta][ic][iRP] = hTmp->ProjectionX(Form("fHistDeltaPtB2_c%d_rp%d",ic,iRP),ibLo,ibUp,"E");
2056 hDeltaPt[iDelta][ic][iRP]->Rebin(2);
2057 if(kNormDelta[iDelta]==2){
2058 if(iRP==0)normRP = hDeltaPt[iDelta][ic][iRP]->Integral("width");
2059 }
2060 else{
2061 normRP = hDeltaPt[iDelta][ic][iRP]->Integral("width");
2062 }
2063 if(normRP) hDeltaPt[iDelta][ic][iRP]->Scale(1./normRP);
2064 }
2065 }
2066 iRP = 0;
2067 for(iRP = 0;iRP<nRP;iRP++){
2068 for(int im = 0;im <nMult;im++){
8f206bf0 2069 TH2F *hTmp = GetTH2PlotB(cBB.Data(),0,1,im,iRP-1); // emb jets
8dc333da 2070 if(!hTmp)Printf("%d %s not found",__LINE__,tmpName.Data());
2071 int ibLo = hTmp->GetYaxis()->FindBin(fMinPt[iDelta]);
2072 int ibUp = hTmp->GetYaxis()->FindBin(fMaxPt[iDelta])-1;
2073 Printf("Line:%d bin %d - %d",__LINE__,ibLo,ibUp);
2074 hDeltaPtMult[iDelta][im][iRP] = hTmp->ProjectionX(Form("fHistDeltaPtB2_%d_mult%d_rp%d",iDelta,im,iRP),ibLo,ibUp,"E");
2075 hDeltaPtMult[iDelta][im][iRP]->Rebin(2);
2076 Float_t fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
2077 if(fScale) hDeltaPtMult[iDelta][im][iRP]->Scale(1./fScale);
2078 }
2079 }
2080 iRP = 0;
2081
2082 iDelta = 4;
1354c5f3 2083 sDelta[iDelta][0] = Form("jet: %2.0f-%2.0f GeV",fMinPt[iDelta],fMaxPt[iDelta]);
5f5e3e77 2084
2085
8dc333da 2086 for(iRP = 0;iRP<nRP;iRP++){
2087 for(int ic = 0;ic < nCen;ic++){
8f206bf0 2088 TH2F *hTmp = GetTH2PlotB(cBB.Data(),1,0,ic,iRP-1); // emb jets
8dc333da 2089 if(!hTmp)Printf("%d %s not found",__LINE__,tmpName.Data());
2090 int ibLo = hTmp->GetYaxis()->FindBin(fMinPt[iDelta]);
2091 int ibUp = hTmp->GetYaxis()->FindBin(fMaxPt[iDelta])-1;
8f206bf0 2092 hDeltaPt[iDelta][ic][iRP] = hTmp->ProjectionX(Form("fHistDeltaPtB1_c%d_rp%d",ic,iRP),ibLo,ibUp,"E");
8dc333da 2093 hDeltaPt[iDelta][ic][iRP]->Rebin(2);
8f206bf0 2094 Float_t fScale = hDeltaPt[iDelta][ic][iRP]->Integral("width");
2095 if(fScale) hDeltaPt[iDelta][ic][iRP]->Scale(1./fScale);
8dc333da 2096 }
2097 }
2098 iRP = 0;
5f5e3e77 2099
8dc333da 2100 for(iRP = 0;iRP<nRP;iRP++){
2101 for(int im = 0;im <nMult;im++){
8f206bf0 2102 TH2F *hTmp = GetTH2PlotB(cBB.Data(),1,1,im,iRP-1); // emb jets
8dc333da 2103 if(!hTmp)Printf("%d %s not found",__LINE__,tmpName.Data());
2104 int ibLo = hTmp->GetYaxis()->FindBin(fMinPt[iDelta]);
2105 int ibUp = hTmp->GetYaxis()->FindBin(fMaxPt[iDelta])-1;
2106 Printf("Line:%d bin %d - %d",__LINE__,ibLo,ibUp);
2107 hDeltaPtMult[iDelta][im][iRP] = hTmp->ProjectionX(Form("fHistDeltaPtB2_%d_mult%d_rp%d",iDelta,im,iRP),ibLo,ibUp,"E");
2108 hDeltaPtMult[iDelta][im][iRP]->Rebin(2);
5f5e3e77 2109 Float_t fScale = hDeltaPtMult[iDelta][im][iRP]->Integral("width");
8dc333da 2110 if(fScale) hDeltaPtMult[iDelta][im][iRP]->Scale(1./fScale);
2111 }
2112 }
8f206bf0 2113
8dc333da 2114 iRP = 0;
2115 // Data from marta
2116 iDelta = 5;
2117 TString fOTFQoff = "~/Dropbox/SharedJets/OnTheFlyOutput/PWG4_JetTasksOutput_110422_p2_OTFQOff_000139107_Total.root";
2118 TFile *fQoff = TFile::Open(fOTFQoff.Data());
1354c5f3 2119 sDelta[iDelta][0] = Form("unquenched jet: %2.0f-%2.0f GeV",fMinPt[iDelta],fMaxPt[iDelta]);
8dc333da 2120 if(fQoff){
2121 for(int ic = 0;ic < nCen;ic++){
2122 tmpName = Form("PWG4_BkgFluct%sCent%dB%d","ANTIKT",ic,iB);
2123 TDirectory *dir = (TDirectory*)fQoff->Get(tmpName.Data());
2124 if(!dir)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
2125 tmpName = Form("taskBkgFluct%sCent%dB%d","ANTIKT",ic,iB);
2126 TList *list = (TList*)dir->Get(tmpName.Data());
2127 if(!list)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
2128 TH3F *h3Tmp = (TH3F*)list->FindObject("fh3PtDeltaPtArea");
2129 if(!h3Tmp)Printf("Line:%d %s not found",__LINE__,"fh3PtDeltaPtArea");
2130 hDeltaPt[iDelta][ic][iRP] = h3Tmp->ProjectionY(Form("hDeltaM%d_%d",ic,iDelta),h3Tmp->GetXaxis()->FindBin(fMinPt[iDelta]),h3Tmp->GetXaxis()->FindBin(fMaxPt[iDelta]),
2131 0,-1,"E");
2132 if(ic<2) hDeltaPt[iDelta][ic][iRP]->Rebin(4);
2133 else if(ic==2) hDeltaPt[iDelta][ic][iRP]->Rebin(2);
2134 Float_t fScale = hDeltaPt[iDelta][ic][iRP]->Integral("width");
2135 if(fScale) hDeltaPt[iDelta][ic][iRP]->Scale(1./fScale);
2136 Printf("Line:%d %p",__LINE__,hDeltaPt[iDelta][ic][iRP]);
2137 }
2138 }
2139 else{
2140 Printf("Could not open %s",fOTFQoff.Data());
2141 }
2142
2143 // marta quenched...
2144 iDelta = 6;
5f5e3e77 2145 TString fOTFQon = "/Users/kleinb/Dropbox/SharedJets/OnTheFlyOutput/LeadingJet/PWG4_JetTasksOutput_110505_p2_OTFQOnL_Merged_v2.root";
8dc333da 2146 TFile *fQon = TFile::Open(fOTFQon.Data());
1354c5f3 2147 sDelta[iDelta][0] = Form("quen. lead. jet: %2.0f-%2.0f GeV",fMinPt[iDelta],fMaxPt[iDelta]);
8dc333da 2148 if(fQon){
2149 for(int ic = 0;ic < nCen;ic++){
2150 tmpName = Form("PWG4_BkgFluct%sCent%dB%d","ANTIKT",ic,iB);
2151 TDirectory *dir = (TDirectory*)fQon->Get(tmpName.Data());
2152 if(!dir)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
2153 tmpName = Form("taskBkgFluct%sCent%dB%d","ANTIKT",ic,iB);
2154 TList *list = (TList*)dir->Get(tmpName.Data());
2155 if(!list)Printf("Line:%d %s not found",__LINE__,tmpName.Data());
2156
2157 TH3F *h3Tmp = (TH3F*)list->FindObject("fh3PtDeltaPtArea");
2158 hDeltaPt[iDelta][ic][iRP] = h3Tmp->ProjectionY(Form("hDeltaM%d_%d",ic,iDelta),h3Tmp->GetXaxis()->FindBin(fMinPt[iDelta]),h3Tmp->GetXaxis()->FindBin(fMaxPt[iDelta]),
2159 0,-1,"E");
2160 if(ic<2) hDeltaPt[iDelta][ic][iRP]->Rebin(4);
2161 else if(ic==2) hDeltaPt[iDelta][ic][iRP]->Rebin(2);
2162 Float_t fScale = hDeltaPt[iDelta][ic][iRP]->Integral("width");
2163 if(fScale)hDeltaPt[iDelta][ic][iRP]->Scale(1./fScale);
2164 }
2165 }
2166 else{
2167 Printf("Could not open %s",fOTFQon.Data());
2168 }
2169
8f206bf0 2170 TString sinputFile = "~/Dropbox/SharedJets/Christian/Files/PWG4_JetTasksOutput_Merge_LHC10h_AOD_110429a.root";
8dc333da 2171 TString sinputDir = "PWG4_spec2_clustersAOD_ANTIKT04_B2_Filter00128_Cut00150_Skip02__0000000000_Class00";
2172 // Plot the real jet spectrum, to be normalized later...
2173 iDelta = 7;
1354c5f3 2174 sDelta[iDelta][0] = "scaled jet spectrum";
8dc333da 2175
2176 if(true){
2177 THnSparseF* fhnJetPtRec = (THnSparseF*)GetObjectFromFile("fhnJetPtRec",sinputFile.Data(),sinputDir.Data());
2178 fhnJetPtRec->GetAxis(0)->SetRange(3,3); // take all jets
2179
2180 // in centrality bins....
2181 for(iRP = 0;iRP<nRP;iRP++){
2182 for(int ic = 0;ic < nCen;ic++){
2183 Int_t icLo = fhnJetPtRec->GetAxis(2)->FindBin(fCentLo[ic]+0.1);
2184 Int_t icUp = fhnJetPtRec->GetAxis(2)->FindBin(fCentUp[ic]-0.1);
2185 fhnJetPtRec->GetAxis(2)->SetRange(icLo,icUp);
2186
2187 // iRP range 0 does rese , this is what we want for all :)
2188 fhnJetPtRec->GetAxis(4)->SetRange(iRP,iRP);
2189
2190 hDeltaPt[iDelta][ic][iRP] = fhnJetPtRec->Projection(1,"E");
2191 hDeltaPt[iDelta][ic][iRP]->SetName(Form("hSpectrumJet_%d_C%d_rp%d",iDelta,ic,iRP));
2192 }
2193 fhnJetPtRec->GetAxis(2)->SetRange(0,0);// reset centrality
2194 for(int im = 0;im < nMult;im++){
2195 Int_t ibTLo = fhnJetPtRec->GetAxis(3)->FindBin(multMin[im]);
2196 Int_t ibTUp = fhnJetPtRec->GetAxis(3)->FindBin(multMax[im]);
2197
2198 fhnJetPtRec->GetAxis(3)->SetRange(ibTLo,ibTUp);
2199 // iRP range 0 does rese , this is what we want for all :)
2200 fhnJetPtRec->GetAxis(4)->SetRange(iRP,iRP);
2201
2202 hDeltaPtMult[iDelta][im][iRP] = fhnJetPtRec->Projection(1,"E");
2203 hDeltaPtMult[iDelta][im][iRP]->SetName(Form("hSpectrumJets_%d_M%d_rp%d",iDelta,im,iRP));
2204 }
2205 }
2206 }
2207 iRP = 0;
2208
2209
8f206bf0 2210 TString sinputFileP = "~/Dropbox/SharedJets/Christian/Files/PWG4_JetTasksOutput_Merge_000139107_AOD_110505a.root";
2211 TString sinputDirP = "PWG4_spec2_clustersAOD_ANTIKT02_B2_Filter00128_Cut00150_Skip00_clustersAOD_ANTIKT02_B0_Filter00128_Cut00150_Skip00_0000000000_Class00";
8dc333da 2212 // fetch the poissonian fluctuations! CKB
2213 iDelta = 8;
2214 sDelta[iDelta][0] = "Poissonian limit";
2215 if(true){
8f206bf0 2216 THnSparseF* fhnTrackPtRec = (THnSparseF*)GetObjectFromFile("fhnTrackPtRec",sinputFileP.Data(),sinputDirP.Data());
8dc333da 2217 fhnTrackPtRec->GetAxis(0)->SetRange(2,2); // take all tracks
5f5e3e77 2218 // iRP range 0 does rese , this is what we wantfor all :)
2219 fhnTrackPtRec->GetAxis(4)->SetRange(iRP,iRP);
2220 THnSparseF* fhnEvent = (THnSparseF*)GetObjectFromFile("fhnEvent",sinputFileP.Data(),sinputDirP.Data());
2221 TH1F *fh1Centrality = (TH1F*)GetObjectFromFile("fh1Centrality",sinputFileP.Data(),sinputDirP.Data());
2222 // Multiplicity
2223 TH2F *fh2MultRec = (TH2F*)GetObjectFromFile("fh2MultRec",sinputFileP.Data(),sinputDirP.Data());
2224 TH1D *h1Mult = (TH1D*)fh2MultRec->ProjectionX("h1Mult");
8dc333da 2225
2226
2227 // in centrality bins....
2228
2229 for(int ic = 0;ic < nCen;ic++){
2230 Int_t icLo = fhnTrackPtRec->GetAxis(2)->FindBin(fCentLo[ic]+0.1);
2231 Int_t icUp = fhnTrackPtRec->GetAxis(2)->FindBin(fCentUp[ic]-0.1);
2232 fhnTrackPtRec->GetAxis(2)->SetRange(icLo,icUp);
2233
2234 TH1D *hMult = fhnEvent->Projection(1,"E");
2235 if(!grSigmaDeltaPtCent[iDelta][iRP]){
2236 grSigmaDeltaPtCent[iDelta][iRP] = new TGraphErrors(nCen);
2237 SetGraphAttributes(grSigmaDeltaPtCent[iDelta][iRP],kMarkerDelta[iDelta],kColorDelta[iDelta]);
2238 }
2239
2240 TH1D *hSpectrumTracks = fhnTrackPtRec->Projection(1,"E");
2241 hSpectrumTracks->SetName(Form("hSpectrumTracks_C%d",ic));
2242
2243 // scale with bin width?
2244
2245 for(int ib = 1;ib<hSpectrumTracks->GetNbinsX();ib++){
2246 Float_t val = hSpectrumTracks->GetBinContent(ib);
2247 hSpectrumTracks->SetBinContent(ib,val/hSpectrumTracks->GetBinWidth(ib));
2248 }
2249
2250 int ibLo = fh1Centrality->FindBin(fCentLo[ic]+0.1);
2251 int ibUp = fh1Centrality->FindBin(fCentUp[ic]-0.1);
2252 Float_t fScale = fh1Centrality->Integral(ibLo,ibUp);
2253 if(fScale>0)hSpectrumTracks->Scale(1./fScale);
5f5e3e77 2254
2255 Double_t sigma = 0;
2256 Double_t sigma_error = 0;
2257 Double_t mult = 0;
2258
2259 Double_t mult2 = h1Mult->GetMean(1);
2260 Double_t mult2_e = h1Mult->GetMean(11); // should be the error on the mult :)
2261
2262 if(iRP==0){// all
2263 sigma = GetPoissonFluctuation(hSpectrumTracks,TMath::Pi()*2.*1.8,
2264 TMath::Pi()*0.4*0.4);
2265 sigma_error = sigma/1000.;
2266 mult = hSpectrumTracks->Integral("width");
2267 }
2268 else{
2269 // uses only 1/3 of acceptance
2270 sigma = GetPoissonFluctuation(hSpectrumTracks,TMath::Pi()*2.*1.8/Float_t(nRP-1),
2271 TMath::Pi()*0.4*0.4);
2272 sigma_error = sigma/1000.;
2273 mult = hSpectrumTracks->Integral("width");
2274 }
2275
8dc333da 2276 Double_t mean = hSpectrumTracks->GetMean();
2277 Double_t rms = hSpectrumTracks->GetRMS();
8dc333da 2278 Double_t cent = (fCentUp[ic]+fCentLo[ic])/2;
2279 Double_t cent_e = (fCentUp[ic]-fCentLo[ic])/2;
2280 grSigmaDeltaPtCent[iDelta][iRP]->SetPoint(ic,cent,sigma);
2281 grSigmaDeltaPtCent[iDelta][iRP]->SetPointError(ic,cent_e,sigma_error);
5f5e3e77 2282
2283
8dc333da 2284 delete hSpectrumTracks;
2285 delete hMult;
2286 }
2287 fhnTrackPtRec->GetAxis(2)->SetRange(0,0);// reset
2288
2289 // in multiplicity bins...
2290 for(iRP = 0;iRP<nRP;iRP++){
2291 for(int im = 0;im < nMult;im++){
2292 Int_t ibELo = fhnEvent->GetAxis(1)->FindBin(multMin[im]);
2293 Int_t ibEUp = fhnEvent->GetAxis(1)->FindBin(multMax[im]);
2294 fhnEvent->GetAxis(1)->SetRange(ibELo,ibEUp);
2295
2296 TH1D *hMult = fhnEvent->Projection(1,"E");
2297
2298 if(!grSigmaDeltaPtMult[iDelta][iRP]){
2299 grSigmaDeltaPtMult[iDelta][iRP] = new TGraphErrors(nCen);
2300 SetGraphAttributes(grSigmaDeltaPtMult[iDelta][iRP],kMarkerDelta[iDelta],kColorDelta[iDelta]);
2301 if(iRP>0)SetGraphAttributes(grSigmaDeltaPtMult[iDelta][iRP],kMarkerRP[iRP],kColorDelta[iDelta]); // change only the marke here
2302 }
2303
2304 // mult range
2305 Int_t ibTLo = fhnTrackPtRec->GetAxis(3)->FindBin(multMin[im]);
2306 Int_t ibTUp = fhnTrackPtRec->GetAxis(3)->FindBin(multMax[im]);
2307 fhnTrackPtRec->GetAxis(3)->SetRange(ibTLo,ibTUp);
2308
2309 // iRP range 0 does rese , this is what we wantfor all :)
2310 fhnTrackPtRec->GetAxis(4)->SetRange(iRP,iRP);
2311 // control
2312 /*
2313 TH2D *h2Tmp = fhnTrackPtRec->Projection(1,4,"E");
2314 h2Tmp->Draw("colz");
2315 c1->Update();
2316 if(getchar()=='q')return 1;
2317 */
2318 TH1D *hSpectrumTracks = fhnTrackPtRec->Projection(1,"E");
2319 hSpectrumTracks->SetName(Form("hSpectrumTracks_M%d_rp%d",im,iRP));
2320
2321 // scale with bin width?
2322 for(int ib = 1;ib<hSpectrumTracks->GetNbinsX();ib++){
2323 Float_t val = hSpectrumTracks->GetBinContent(ib);
2324 hSpectrumTracks->SetBinContent(ib,val/hSpectrumTracks->GetBinWidth(ib));
2325 }
2326
2327 int ibLo = h1Mult->FindBin(multMin[im]);
2328 int ibUp = h1Mult->FindBin(multMax[im]);
2329 h1Mult->GetXaxis()->SetRange(ibLo,ibUp);
2330 Float_t fScale = h1Mult->Integral(ibLo,ibUp);
2331 if(fScale>0)hSpectrumTracks->Scale(1./fScale);
2332
2333 Double_t sigma = 0;
2334 Double_t sigma_error = 0;
2335 Double_t mult = 0;
2336
2337 Double_t mult2 = h1Mult->GetMean(1);
2338 Double_t mult2_e = h1Mult->GetMean(11); // should be the error on the mult :)
2339
2340 if(iRP==0){// all
2341 sigma = GetPoissonFluctuation(hSpectrumTracks,TMath::Pi()*2.*1.8,
2342 TMath::Pi()*0.4*0.4);
2343 sigma_error = sigma/1000.;
2344 mult = hSpectrumTracks->Integral("width");
2345 }
2346 else{
2347 // uses only 1/3 of acceptance
2348 sigma = GetPoissonFluctuation(hSpectrumTracks,TMath::Pi()*2.*1.8/Float_t(nRP-1),
2349 TMath::Pi()*0.4*0.4);
2350 sigma_error = sigma/1000.;
2351 mult = hSpectrumTracks->Integral("width");
2352 }
2353 Printf("mult %4.3f mult2 %4.3f",mult,mult2);
2354
2355 grSigmaDeltaPtMult[iDelta][iRP]->SetPoint(im,mult2,sigma);
2356 grSigmaDeltaPtMult[iDelta][iRP]->SetPointError(im,mult2_e,sigma_error);
5f5e3e77 2357
2358
2359 // correct for eff
2360 Int_t oldDelta = iDelta;
2361 iDelta = 11;
2362
2363 sDelta[iDelta][0] = sDelta[oldDelta][0] + " eff. corrected";
2364
2365 if(!grSigmaDeltaPtMult[iDelta][iRP]){
2366 grSigmaDeltaPtMult[iDelta][iRP] = new TGraphErrors(nCen);
2367 SetGraphAttributes(grSigmaDeltaPtMult[iDelta][iRP],kMarkerDelta[iDelta],kColorDelta[iDelta]);
2368 if(iRP>0)SetGraphAttributes(grSigmaDeltaPtMult[iDelta][iRP],kMarkerRP[iRP],kColorDelta[iDelta]); // change only the marke here
2369 }
2370
2371 TH1* hCorr = CorrectForEff(hSpectrumTracks);
2372
2373 if(iRP==0){// all
2374 sigma = GetPoissonFluctuation(hCorr,TMath::Pi()*2.*1.8,
2375 TMath::Pi()*0.4*0.4);
2376 sigma_error = sigma/1000.;
2377 mult = hCorr->Integral("width");
2378 }
2379 else{
2380 // uses only 1/3 of acceptance
2381 sigma = GetPoissonFluctuation(hCorr,TMath::Pi()*2.*1.8/Float_t(nRP-1),
2382 TMath::Pi()*0.4*0.4);
2383 sigma_error = sigma/1000.;
2384 mult = hCorr->Integral("width");
2385 }
2386 Printf("mult %4.3f mult2 %4.3f",mult,mult2);
2387
2388 grSigmaDeltaPtMult[iDelta][iRP]->SetPoint(im,mult2,sigma);
2389 grSigmaDeltaPtMult[iDelta][iRP]->SetPointError(im,mult2_e,sigma_error);
2390
2391
2392 iDelta = oldDelta;
8dc333da 2393 delete hSpectrumTracks;
2394 delete hMult;
2395 }
2396 }
2397 // in multiplicity bins...
2398
2399
2400 }
2401
8f206bf0 2402 for(iDelta = 0;iDelta < kDeltaTypes;iDelta++){
1354c5f3 2403 sDelta[iDelta][0] += (iPlotType&(1<<2))?" (all)":"";
2404 sDelta[iDelta][1] = "in plane";
2405 sDelta[iDelta][2] = "in between";
2406 sDelta[iDelta][3] = "out of plane";
8dc333da 2407 }
2408
2409 // THIS IS FOR THE ALICE LABELS AND THE WORK IN PROGRESS PRELIMINARY ETC.
2410 TLatex *txt2 = new TLatex();
2411 txt2->SetTextFont(gStyle->GetTextFont());
2412 txt2->SetTextSize(gStyle->GetTextSize()*0.7);
2413 txt2->SetTextAlign(22);
2414 txt2->SetTextColor(kRed);
1354c5f3 2415 TLatex* txtHead = new TLatex();
2416 txtHead->SetTextSize(0.04);
2417 txtHead->SetNDC();
2418 txtHead->SetTextAlign(23);
8dc333da 2419
2420 if(iPlotType&(1<<0)){
2421 for(int ic = 0;ic < nCen;ic++){
2422 TLegend *leg1 = new TLegend(0.2,0.7,0.3,0.98);
8dc333da 2423 leg1->SetFillColor(0);
2424 leg1->SetTextFont(gStyle->GetTextFont());
2425 leg1->SetTextSize(gStyle->GetTextSize()*0.6);
2426 leg1->SetBorderSize(0);
2427 hFrame->DrawCopy();
2428
1354c5f3 2429 Int_t iLegCount = -1;
8dc333da 2430 for(iRP = 0;iRP<((iPlotType&(1<<2))?nRP:1);iRP++){
2431 Double_t mean = 0;
2432 Double_t sigma = 0;
2433 Double_t sigma_error = 0;
2434 Double_t mean_error = 0;
2435 TF1* tmpGaus = 0;
2436 for(iDelta = 0;iDelta<kDeltaTypes;iDelta++){
5f5e3e77 2437
8dc333da 2438 if(!grMeanDeltaPtCent[iDelta][iRP]&&kFitDelta[iDelta]){
2439 grMeanDeltaPtCent[iDelta][iRP] = new TGraphErrors(nCen);
2440 SetGraphAttributes(grMeanDeltaPtCent[iDelta][iRP],kMarkerDelta[iDelta],kColorDelta[iDelta]);
2441 if(iRP>0)SetGraphAttributes(grMeanDeltaPtCent[iDelta][iRP],kMarkerRP[iRP],kColorDelta[iRP]);
2442 }
2443 if(!grSigmaDeltaPtCent[iDelta][iRP]&&kFitDelta[iDelta]){
2444 grSigmaDeltaPtCent[iDelta][iRP] = new TGraphErrors(nCen);
2445 SetGraphAttributes(grSigmaDeltaPtCent[iDelta][iRP],kMarkerDelta[iDelta],kColorDelta[iDelta]);
2446 if(iRP>0)SetGraphAttributes(grSigmaDeltaPtCent[iDelta][iRP],kMarkerRP[iRP],kColorDelta[iRP]);
2447 }
2448
5f5e3e77 2449 if(!hDeltaPt[iDelta][ic][iRP]){
2450 Printf("%d:%d:%d not found",iDelta,ic,iRP);
2451 continue;
2452 }
2453 SetHistoAttributes(hDeltaPt[iDelta][ic][iRP],kMarkerDelta[iDelta],kColorDelta[iDelta]);
2454 if(iRP!=0)SetHistoAttributes(hDeltaPt[iDelta][ic][iRP],kMarkerRP[iRP],kColorRP[iRP]);
2455 if(!(iPlotFlag&(1<<iDelta)))continue;
2456
2457
8dc333da 2458 // take the first iDelta Pt as anchor for norm, worst case does nothin (scale with 1
2459 if(kNormValue[ic][iRP]<=0){
2460 Int_t ib = hDeltaPt[iDelta][ic][iRP]->FindBin(kNormPt);
2461 kNormValue[ic][iRP] = hDeltaPt[iDelta][ic][iRP]->GetBinContent(ib);
2462 }
2463 if(iDelta==kNormType){
2464 Int_t ib = hDeltaPt[iDelta][ic][iRP]->FindBin(kNormPt);
2465 Float_t val1 = hDeltaPt[iDelta][ic][iRP]->GetBinContent(ib);
2466 if(val1){
2467 hDeltaPt[iDelta][ic][iRP]->Scale(kNormValue[ic][iRP]/val1);
2468 }
2469 }
1354c5f3 2470 iLegCount++;
2471 hDeltaPt[iDelta][ic][iRP]->DrawCopy("psame");
2472 // leg1->AddEntry(hDeltaPt[iDelta][ic][iRP],sDelta[iDelta][iRP].Data(),"P");
2473 c1->Update();
2474 draw_legend_m(0.01,0.95,(TAttMarker*)hDeltaPt[iDelta][ic][iRP],sDelta[iDelta][iRP].Data(),iLegCount,1,0.026);
8dc333da 2475
2476
8dc333da 2477 if(kFitDelta[iDelta]){
2478 tmpGaus = FitLHSgaus(hDeltaPt[iDelta][ic][iRP]);
2479 mean = tmpGaus->GetParameter(1);
2480 sigma = tmpGaus->GetParameter(2);
2481 mean_error = tmpGaus->GetParError(1);
2482 sigma_error = tmpGaus->GetParError(2);
2483
2484 tmpGaus->SetRange(-40,40);
8f206bf0 2485 tmpGaus->SetLineColor(hDeltaPt[iDelta][ic][iRP]->GetLineColor());
1354c5f3 2486 tmpGaus->SetLineStyle(kLineDelta[iDelta]);
8dc333da 2487 tmpGaus->Draw("same");
1354c5f3 2488 // leg1->AddEntry(tmpGaus,Form("LHS Gauss fit: #mu = %2.2f, #sigma = %2.2f",mean,sigma),"L");
2489 draw_legend_l(0.45,0.95,tmpGaus,Form("LHS fit: #mu = %2.2f, #sigma = %2.2f",mean,sigma),iLegCount,1,0.026);
8dc333da 2490 Double_t cent = (fCentUp[ic]+fCentLo[ic])/2;
2491 Double_t cent_e = (fCentUp[ic]-fCentLo[ic])/2;
2492 Printf("cent %3.3f +- %3.3f",cent,cent_e);
2493 grMeanDeltaPtCent[iDelta][iRP]->SetPoint(ic,cent,mean);
2494 grMeanDeltaPtCent[iDelta][iRP]->SetPointError(ic,cent_e,mean_error);
2495
2496 grSigmaDeltaPtCent[iDelta][iRP]->SetPoint(ic,cent,sigma);
2497 grSigmaDeltaPtCent[iDelta][iRP]->SetPointError(ic,cent_e,sigma_error);
2498 }
2499 } //deltatypes
2500 }// RP
1354c5f3 2501 // leg1->Draw();
2502 txtHead->DrawLatex(0.5,0.99,Form("LHC2010 Pb-Pb %2.0f-%2.0f%% R = 0.4 (B%d)",fCentLo[ic],fCentUp[ic],iB));
2503 DrawALICE(c1,iALICEType,0.7,0.72);
8dc333da 2504 c1->Update();
2505 c1->SaveAs(Form("%sdeltaPt_pT%s%s_B%d_cen%02d_%03d.%s",
2506 picPrefix.Data(),
2507 Form("%03.0f_%03.0f",inPtLo,inPtUp),
2508 (iPlotType&(1<<2))?"_rp":"",
2509 iB,ic,iPlotFlag,
2510 picSuffix.Data()));
2511 if(!gROOT->IsBatch()){
2512 if(getchar()=='q')return 1;
2513 }
2514 }
2515
2516 // Draw the trending plots vs cent...
2517 c1->SetLogy(0);
1354c5f3 2518 TH2F* hFrameMeanCent = new TH2F("hMeanCent",";centrality (%);#mu of LHS Gauss (GeV/c)",1,0.,100.,1,-10.,15.);
8dc333da 2519 hFrameMeanCent->SetLabelSize(hFrame->GetLabelSize("Y")*0.9,"Y");
2520 hFrameMeanCent->SetLabelSize(hFrame->GetLabelSize("X")*0.9,"X");
2521 hFrameMeanCent->SetTitleSize(hFrame->GetTitleSize("Y")*0.9,"Y");
2522 hFrameMeanCent->SetTitleSize(hFrame->GetTitleSize("X")*0.9,"X");
2523 hFrameMeanCent->SetTitleOffset(1.1,"Y");
2524 hFrameMeanCent->SetTitleOffset(1.1,"X");
2525 hFrameMeanCent->SetStats(kFALSE);
2526 hFrameMeanCent->DrawCopy();
2527
2528 TLegend *legSM = new TLegend(0.2,0.7,0.3,0.98);
1354c5f3 2529 // legSM->SetHeader(Form("Pb-Pb R = 0.4 (B%d)",iB));
8dc333da 2530 legSM->SetFillColor(0);
2531 legSM->SetTextFont(gStyle->GetTextFont());
2532 legSM->SetTextSize(gStyle->GetTextSize()*0.6);
2533 legSM->SetBorderSize(0);
2534
1354c5f3 2535 Int_t iLegCount = 0;
8dc333da 2536 for(iRP = 0;iRP<((iPlotType&(1<<2))?nRP:1);iRP++){
2537 for(iDelta = 0;iDelta <kDeltaTypes;iDelta++){
2538 if(!(iPlotFlag&(1<<iDelta)))continue;
2539 if(!grMeanDeltaPtCent[iDelta][iRP])continue;
1354c5f3 2540 c1->Update();
2541 if(iDelta==8||iDelta==8){
2542 }
2543 else{
2544 grMeanDeltaPtCent[iDelta][iRP]->Draw("psame");
2545 draw_legend_m(0.01,0.95,grMeanDeltaPtCent[iDelta][iRP],sDelta[iDelta][iRP].Data(),iLegCount,1,0.026);
2546 // legSM->AddEntry(grMeanDeltaPtCent[iDelta][iRP],sDelta[iDelta][iRP].Data(),"P");
2547 iLegCount++;
2548 }
8dc333da 2549 }
2550 }
1354c5f3 2551 // legSM->Draw();
2552 txtHead->DrawLatex(0.5,0.99,Form("LHC2010 Pb-Pb #sqrt{s} = 2.76 TeV"));
2553 DrawALICE(c1,iALICEType,0.7,0.8);
8dc333da 2554 c1->Update();
2555 c1->SaveAs(Form("%sMeanVsCent_pT%s%s_B%d_%03d.%s",
2556 picPrefix.Data(),
2557 Form("%03.0f_%03.0f",inPtLo,inPtUp),
2558 (iPlotType&(1<<2))?"_rp":"",
2559 iB,iPlotFlag,picSuffix.Data()));
2560 if(!gROOT->IsBatch()){
2561 if(getchar()=='q')return 1;
2562 }
1354c5f3 2563 TH2F* hFrameSigmaCent = new TH2F("hSigmaCent",";centrality (%);#sigma of LHS Gauss (GeV/c)",1,0.,100.,10,0.,18.);
8dc333da 2564 hFrameSigmaCent->SetLabelSize(hFrame->GetLabelSize("Y")*0.9,"Y");
2565 hFrameSigmaCent->SetLabelSize(hFrame->GetLabelSize("X")*0.9,"X");
2566 hFrameSigmaCent->SetTitleSize(hFrame->GetTitleSize("Y")*0.9,"Y");
2567 hFrameSigmaCent->SetTitleSize(hFrame->GetTitleSize("X")*0.9,"X");
2568 hFrameSigmaCent->SetTitleOffset(1.1,"Y");
2569 hFrameSigmaCent->SetTitleOffset(1.1,"X");
2570 hFrameSigmaCent->SetStats(kFALSE);
2571 hFrameSigmaCent->DrawCopy();
2572
1354c5f3 2573 /*
8dc333da 2574 TLegend *legSC = new TLegend(0.2,0.7,0.3,0.98);
1354c5f3 2575 legSC->SetHeader(Form("Pb-Pb R = 0.4 (B%d)",iB));
8dc333da 2576 legSC->SetFillColor(0);
2577 legSC->SetTextFont(gStyle->GetTextFont());
2578 legSC->SetTextSize(gStyle->GetTextSize()*0.6);
2579 legSC->SetBorderSize(0);
1354c5f3 2580 */
2581
2582 iLegCount = 0;
2583
8dc333da 2584 for(iRP = 0;iRP<((iPlotType&(1<<2))?nRP:1);iRP++){
2585 for(iDelta = 0;iDelta <kDeltaTypes;iDelta++){
2586 if(!(iPlotFlag&(1<<iDelta)))continue;
2587 if(!grSigmaDeltaPtCent[iDelta][iRP])continue;
2588 grSigmaDeltaPtCent[iDelta][iRP]->Draw("psame");
1354c5f3 2589 // legSC->AddEntry(grSigmaDeltaPtCent[iDelta][iRP],sDelta[iDelta][iRP].Data(),"P");
2590 c1->Update();
2591 draw_legend_m(0.01,0.95,grSigmaDeltaPtCent[iDelta][iRP],sDelta[iDelta][iRP].Data(),iLegCount,1,0.026);
2592 iLegCount++;
8dc333da 2593 }
2594 }
1354c5f3 2595 // legSC->Draw();
2596 txtHead->DrawLatex(0.5,0.99,Form("LHC2010 Pb-Pb #sqrt{s} = 2.76 TeV"));
2597 DrawALICE(c1,iALICEType,0.7,0.8);
8dc333da 2598 c1->Update();
8dc333da 2599 c1->SaveAs(Form("%sSigmaVsCent_pT%s%s_B%d_%03d.%s",
2600 picPrefix.Data(),
2601 Form("%03.0f_%03.0f",inPtLo,inPtUp),
2602 (iPlotType&(1<<2))?"_rp":"",
2603 iB,iPlotFlag,picSuffix.Data()));
2604
2605 if(!gROOT->IsBatch()){
2606 if(getchar()=='q')return 1;
2607 }
2608 }
2609
2610 // plot the trends vs mult
2611 if(iPlotType&(1<<1)){
2612 for(int im = 0;im < nMult;im++){
2613 TLegend *leg1 = new TLegend(0.2,0.7,0.3,0.98);
1354c5f3 2614
8dc333da 2615 leg1->SetFillColor(0);
2616 leg1->SetTextFont(gStyle->GetTextFont());
2617 leg1->SetTextSize(gStyle->GetTextSize()*0.6);
2618 leg1->SetBorderSize(0);
1354c5f3 2619 c1->SetLogy();
8dc333da 2620 hFrame->DrawCopy();
2621
1354c5f3 2622 Int_t iLegCount = -1;
8dc333da 2623 for(iRP = 0;iRP<((iPlotType&(1<<2))?nRP:1);iRP++){
2624 Double_t mean = 0;
2625 Double_t sigma = 0;
2626 Double_t sigma_error = 0;
2627 Double_t mean_error = 0;
2628 TF1* tmpGaus = 0;
2629 for(iDelta = 0;iDelta<kDeltaTypes;iDelta++){
8dc333da 2630
2631 if(!grMeanDeltaPtMult[iDelta][iRP]&&kFitDelta[iDelta]){
2632 grMeanDeltaPtMult[iDelta][iRP] = new TGraphErrors(nMult);
8dc333da 2633 }
2634 if(!grSigmaDeltaPtMult[iDelta][iRP]&&kFitDelta[iDelta]){
2635 grSigmaDeltaPtMult[iDelta][iRP] = new TGraphErrors(nMult);
8dc333da 2636
2637 }
2638
5f5e3e77 2639 SetGraphAttributes(grMeanDeltaPtMult[iDelta][iRP],kMarkerDelta[iDelta],kColorDelta[iDelta]);
2640 if(iRP>0)SetGraphAttributes(grMeanDeltaPtMult[iDelta][iRP],kMarkerRP[iRP],kColorRP[iRP]);
2641
2642 SetGraphAttributes(grSigmaDeltaPtMult[iDelta][iRP],kMarkerDelta[iDelta],kColorDelta[iDelta]);
2643 if(iRP>0)SetGraphAttributes(grSigmaDeltaPtMult[iDelta][iRP],kMarkerRP[iRP],kColorRP[iRP]);
2644
2645 if(!hDeltaPtMult[iDelta][im][iRP]){
2646 Printf("%d:%d:%d not found",iDelta,im,iRP);
2647 continue;
2648 }
2649 SetHistoAttributes(hDeltaPtMult[iDelta][im][iRP],kMarkerDelta[iDelta],kColorDelta[iDelta]);
2650 if(iRP>0)SetHistoAttributes(hDeltaPtMult[iDelta][im][iRP],kMarkerRP[iRP],kColorRP[iRP]);
2651 if(!(iPlotFlag&(1<<iDelta)))continue;
2652
2653
8dc333da 2654 // take the first iDelta Pt as anchor for norm, worst case does nothin (scale with 1
2655 if(kNormValueMult[im][iRP]<=0){
2656 Int_t ib = hDeltaPtMult[iDelta][im][iRP]->FindBin(kNormPt);
2657 kNormValueMult[im][iRP] = hDeltaPtMult[iDelta][im][iRP]->GetBinContent(ib);
2658 }
2659
2660
2661 if(iDelta==kNormType){
2662 Int_t ib = hDeltaPtMult[iDelta][im][iRP]->FindBin(kNormPt);
2663 Float_t val1 = hDeltaPtMult[iDelta][im][iRP]->GetBinContent(ib);
2664 if(val1!=0){
2665 hDeltaPtMult[iDelta][im][iRP]->Scale(kNormValueMult[im][iRP]/val1);
2666 }
2667 }
2668
2669 hDeltaPtMult[iDelta][im][iRP]->DrawCopy("psame");
2670 leg1->AddEntry(hDeltaPtMult[iDelta][im][iRP],sDelta[iDelta][iRP].Data(),"P");
1354c5f3 2671 c1->Update();
2672 iLegCount++;
2673 draw_legend_m(0.01,0.95,hDeltaPtMult[iDelta][im][iRP],sDelta[iDelta][iRP].Data(),iLegCount,1,0.026);
8dc333da 2674 if(kFitDelta[iDelta]){
2675 tmpGaus = FitLHSgaus(hDeltaPtMult[iDelta][im][iRP]);
2676 mean = tmpGaus->GetParameter(1);
2677 sigma = tmpGaus->GetParameter(2);
2678 mean_error = tmpGaus->GetParError(1);
2679 sigma_error = tmpGaus->GetParError(2);
2680
1354c5f3 2681 /*
2682 if(iDelta==8){
2683 mean= hDeltaPtMult[iDelta][im][iRP]->GetMean();
2684 sigma= hDeltaPtMult[iDelta][im][iRP]->GetRMS();
2685 }
2686 */
8dc333da 2687 tmpGaus->SetRange(-40,40);
2688 tmpGaus->SetLineColor(hDeltaPtMult[iDelta][im][iRP]->GetLineColor());
2689 tmpGaus->Draw("same");
1354c5f3 2690 leg1->AddEntry(tmpGaus,Form("LHS Gauss fit: #mu = %2.2f, #sigma = %2.2f",mean,sigma),"L");
2691 draw_legend_l(0.45,0.95,tmpGaus,sDelta[iDelta][iRP].Data(),iLegCount,1,0.026);
8dc333da 2692 // to be replaced by actual mean of the mutliplicity
2693 Double_t mult = (multMin[im]+multMax[im])/2;
2694 Double_t mult_e = (multMax[im]-multMin[im])/2;
2695 Printf("mult %3.3f +- %3.3f",mult,mult_e);
2696 grMeanDeltaPtMult[iDelta][iRP]->SetPoint(im,mult,mean);
2697 grMeanDeltaPtMult[iDelta][iRP]->SetPointError(im,mult_e,mean_error);
2698
8f206bf0 2699
2700 grSigmaDeltaPtMult[iDelta][iRP]->SetPoint(im,mult,sigma);
2701 grSigmaDeltaPtMult[iDelta][iRP]->SetPointError(im,mult_e,sigma_error);
2702
8dc333da 2703 }
2704 }
2705 }
1354c5f3 2706 // leg1->Draw();
2707 txtHead->DrawLatex(0.5,0.99,Form("Pb-Pb mult. %03d-%03d R = 0.4 (B%d)",multMin[im],multMax[im],iB));
8dc333da 2708 txt2->SetNDC();
1354c5f3 2709
8dc333da 2710 // txt2->DrawLatex(0.5,0.97,"ALICE Performance 01/03/2011");
1354c5f3 2711 DrawALICE(c1,iALICEType,0.7,0.72);
8dc333da 2712 c1->Update();
2713 c1->SaveAs(Form("%sdeltaPt_pT%s%s_B%d_mult%02d_%03d.%s",
2714 picPrefix.Data(),
2715 Form("%03.0f_%03.0f",inPtLo,inPtUp),
2716 (iPlotType&(1<<2))?"_rp":"",
2717 iB,im,iPlotFlag,
2718 picSuffix.Data()));
2719 if(!gROOT->IsBatch()){
2720 if(getchar()=='q')return 1;
2721 }
2722 }
2723 // Draw the trending plots vs cent...
2724 c1->SetLogy(0);
1354c5f3 2725 TH2F* hFrameMeanMult = new TH2F("hMeanMult",";raw #tracks;#mu of LHS Gauss (GeV/c)",500,0.,3200.,12,-10.,15.);
8dc333da 2726 hFrameMeanMult->SetLabelSize(hFrame->GetLabelSize("Y")*0.9,"Y");
2727 hFrameMeanMult->SetLabelSize(hFrame->GetLabelSize("X")*0.9,"X");
2728 hFrameMeanMult->SetTitleSize(hFrame->GetTitleSize("Y")*0.9,"Y");
2729 hFrameMeanMult->SetTitleSize(hFrame->GetTitleSize("X")*0.9,"X");
2730 hFrameMeanMult->SetTitleOffset(1.1,"Y");
2731 hFrameMeanMult->SetTitleOffset(1.1,"X");
2732 hFrameMeanMult->SetStats(kFALSE);
2733 hFrameMeanMult->DrawCopy();
2734
2735 TLegend *legMM = new TLegend(0.2,0.7,0.3,0.98);
1354c5f3 2736 legMM->SetHeader(Form("Pb-Pb R = 0.4 (B%d)",iB));
8dc333da 2737 legMM->SetFillColor(0);
2738 legMM->SetTextFont(gStyle->GetTextFont());
2739 legMM->SetTextSize(gStyle->GetTextSize()*0.6);
2740 legMM->SetBorderSize(0);
2741
1354c5f3 2742 iLegCount = -1;
2743 c1->Update();
8dc333da 2744 for(iRP = 0;iRP<((iPlotType&(1<<2))?nRP:1);iRP++){
2745 for(iDelta = 0;iDelta <kDeltaTypes;iDelta++){
2746 if(!(iPlotFlag&(1<<iDelta)))continue;
2747 if(!grMeanDeltaPtMult[iDelta][iRP])continue;
8f206bf0 2748
5f5e3e77 2749 if(iDelta==8||iDelta==11){
1354c5f3 2750 // iLegCount++;
2751 // grMeanDeltaPtMult[iDelta][iRP]->Draw("Csame");
2752 // legMM->AddEntry(grSigmaDeltaPtMult[iDelta][iRP],sDelta[iDelta][iRP].Data(),"L");
2753 // draw_legend_l(0.01,0.95,grMeanMultPtCent[iDelta][iRP],sDelta[iDelta][iRP].Data(),iLegCount,1,0.026);
8f206bf0 2754 }
2755 else{
1354c5f3 2756 iLegCount++;
8f206bf0 2757 grMeanDeltaPtMult[iDelta][iRP]->Draw("psame");
1354c5f3 2758 // legMM->AddEntry(grSigmaDeltaPtMult[iDelta][iRP],sDelta[iDelta][iRP].Data(),"P");
2759 draw_legend_m(0.01,0.95,grMeanDeltaPtMult[iDelta][iRP],sDelta[iDelta][iRP].Data(),iLegCount,1,0.026);
8f206bf0 2760 }
2761
8dc333da 2762 }
2763 }
1354c5f3 2764 // legMM->Draw();
2765 txtHead->DrawLatex(0.5,0.99,Form("LHC2010 Pb-Pb #sqrt{s} = 2.76 TeV"));
2766 DrawALICE(c1,iALICEType,0.7,0.8);
2767
8dc333da 2768 c1->Update();
8dc333da 2769 c1->SaveAs(Form("%sMeanVsMult_pT%s%s_B%d_%03d.%s",
2770 picPrefix.Data(),
2771 Form("%03.0f_%03.0f",inPtLo,inPtUp),
2772 (iPlotType&(1<<2))?"_rp":"",
2773 iB,iPlotFlag,picSuffix.Data()));
2774 if(!gROOT->IsBatch()){
2775 if(getchar()=='q')return 1;
2776 }
1354c5f3 2777 TH2F* hFrameSigmaMult = new TH2F("hSigmaMult",";raw #tracks;#sigma of LHS Gauss (GeV/c)",500,0.,3200.,10,0.,18.);
8dc333da 2778 hFrameSigmaMult->SetLabelSize(hFrame->GetLabelSize("Y")*0.9,"Y");
2779 hFrameSigmaMult->SetLabelSize(hFrame->GetLabelSize("X")*0.9,"X");
2780 hFrameSigmaMult->SetTitleSize(hFrame->GetTitleSize("Y")*0.9,"Y");
2781 hFrameSigmaMult->SetTitleSize(hFrame->GetTitleSize("X")*0.9,"X");
2782 hFrameSigmaMult->SetTitleOffset(1.1,"Y");
2783 hFrameSigmaMult->SetTitleOffset(1.1,"X");
2784 hFrameSigmaMult->SetStats(kFALSE);
2785 hFrameSigmaMult->DrawCopy();
2786
2787 TLegend *legSM = new TLegend(0.2,0.7,0.3,0.98);
1354c5f3 2788 legSM->SetHeader(Form("Pb-Pb R = 0.4 (B%d)",iB));
8dc333da 2789 legSM->SetFillColor(0);
2790 legSM->SetTextFont(gStyle->GetTextFont());
2791 legSM->SetTextSize(gStyle->GetTextSize()*0.6);
2792 legSM->SetBorderSize(0);
1354c5f3 2793
2794 iLegCount=-1;
2795 c1->Update();
5f5e3e77 2796 for(iRP = ((iPlotType&(1<<2))?1:0);iRP<((iPlotType&(1<<2))?nRP:1);iRP++){
8dc333da 2797 for(iDelta = 0;iDelta <kDeltaTypes;iDelta++){
2798 if(!(iPlotFlag&(1<<iDelta)))continue;
2799 if(!grSigmaDeltaPtMult[iDelta][iRP])continue;
1354c5f3 2800 Printf("%d leg coutn",iLegCount);
5f5e3e77 2801 if(iDelta==8||iDelta==11){
1354c5f3 2802 iLegCount++;
5f5e3e77 2803 grSigmaDeltaPtMult[iDelta][iRP]->Draw("Csame");
1354c5f3 2804 TString cTmp1 = sDelta[iDelta][0];
2805 cTmp1.ReplaceAll(" (all)"," ");
2806 sDelta[iDelta][iRP] = cTmp1 + sDelta[iDelta][iRP];
2807 draw_legend_l(0.01,0.95,grSigmaDeltaPtMult[iDelta][iRP],sDelta[iDelta][iRP].Data(),iLegCount,1,0.026);
5f5e3e77 2808 }
2809 else{
1354c5f3 2810 iLegCount++;
2811 TString cTmp1 = sDelta[iDelta][0];
2812 cTmp1.ReplaceAll(" (all)"," ");
2813 sDelta[iDelta][iRP] = cTmp1 + sDelta[iDelta][iRP];
5f5e3e77 2814 grSigmaDeltaPtMult[iDelta][iRP]->Draw("psame");
1354c5f3 2815 draw_legend_m(0.01,0.95,grSigmaDeltaPtMult[iDelta][iRP],sDelta[iDelta][iRP].Data(),iLegCount,1,0.026);
5f5e3e77 2816 }
1354c5f3 2817
8dc333da 2818 }
2819 }
1354c5f3 2820 // legSM->Draw();
2821 txtHead->DrawLatex(0.5,0.99,Form("LHC2010 Pb-Pb #sqrt{s} = 2.76 TeV"));
2822
2823 DrawALICE(c1,iALICEType,0.7,0.8);
8dc333da 2824 c1->Update();
2825 c1->SaveAs(Form("%sSigmaVsMult_pT%s%s_B%d_%03d.%s",
2826 picPrefix.Data(),
2827 Form("%03.0f_%03.0f",inPtLo,inPtUp),
2828 (iPlotType&(1<<2))?"_rp":"",
2829 iB,iPlotFlag,picSuffix.Data()));
2830 if(!gROOT->IsBatch()){
2831 if(getchar()=='q')return 1;
2832 }
2833
2834 }// trends vs. mult
2835
2836 // plot the trend with RP and vs. mult
2837 CloseFiles();
2838 return 0;
2839
2840
2841}
2842
2843
1354c5f3 2844void ScaleH1(TH1* h,char* cX,char* cY,Float_t fScale,Bool_t bWidth){
e2081c22 2845 if(!h)return;
2846 if(!fScale)return;
2847
2848 h->Scale(fScale);
2849 if(bWidth){
2850 for(int ib = 1;ib <= h->GetNbinsX();ib++){
2851 Float_t val = h->GetBinContent(ib);
2852 Float_t err = h->GetBinError(ib);
2853 Float_t w = h->GetBinWidth(ib);
2854 h->SetBinContent(ib,val/w);
2855 h->SetBinError(ib,err/w);
2856 // Printf("width %f",w);
2857 }
2858 }
2859 h->SetXTitle(cX);
2860 h->SetYTitle(cY);
2861
2862}
2863
2864
2865void ScaleH2(TH2* h,char* cX,char* cY,char* cZ,Float_t fScale,Bool_t bWidth){
2866 if(!h)return;
2867 if(!fScale)return;
2868
2869 h->Scale(fScale);
2870 if(bWidth){
2871 for(int ibx = 1;ibx <= h->GetNbinsX();ibx++){
2872 for(int iby = 1;iby <= h->GetNbinsY();iby++){
2873 Float_t val = h->GetBinContent(ibx,iby);
2874 Float_t err = h->GetBinError(ibx,iby);
2875 Float_t wx = h->GetXaxis()->GetBinWidth(ibx);
2876 Float_t wy = h->GetYaxis()->GetBinWidth(iby);
2877 h->SetBinContent(ibx,iby,val/(wx*wy));
2878 h->SetBinError(ibx,iby,err/(wx*wy));
2879 // Printf("width %f",w);
2880 }
2881 }
2882 }
2883 h->SetXTitle(cX);
2884 h->SetYTitle(cY);
2885 h->SetZTitle(cZ);
2886
2887}
2888
2889void SetStyleH1(TH1 *h,Int_t color,Int_t fillStyle,Int_t markerStyle){
2890 if(color){
2891 h->SetLineColor(color);
2892 h->SetLineWidth(3);
2893 h->SetMarkerColor(color);
2894 h->SetFillColor(color);
2895 }
2896 if(fillStyle)h->SetFillStyle(fillStyle);
2897 if(markerStyle)h->SetMarkerStyle(markerStyle);
2898}
2899
2900void HistsFromSingleJets(char* cMask,TList *list,Int_t nJets,TH1F** h,Int_t iRebin,int iFirst){
2901
2902 for(int ij = iFirst;ij < nJets;ij++){
2903 h[ij] = (TH1F*)list->FindObject(Form(cMask,ij));
2904
2905 if(h[ij]){
2906 if(iRebin>1){
2907 h[ij] = (TH1F*)h[ij]->Clone();
2908 h[ij]->SetDirectory(gROOT);
2909 h[ij]->Rebin(iRebin);
2910 }
2911 if(ij==iFirst){
2912 h[nJets] = (TH1F*)h[ij]->Clone(Form(cMask,nJets));
2913 }
2914 else{
2915 h[nJets]->Add(h[ij]);
2916 }
2917 }
2918 else{
2919 Printf("%s not found",Form(cMask,ij));
2920 }
2921 }
2922
2923}
2924
2925
2926void HistsFromSingleJets(char* cMask,TList *list,Int_t nJets,TH2F** h,Int_t iRebin,int iFirst){
2927
2928 for(int ij = iFirst;ij < nJets;ij++){
2929 h[ij] = (TH2F*)list->FindObject(Form(cMask,ij));
2930 if(h[ij]){
2931 if(iRebin>1){
2932 h[ij] = (TH2F*)h[ij]->Clone();
2933 h[ij]->SetDirectory(gROOT);
2934 h[ij]->Rebin(iRebin);
2935 }
2936 if(ij==iFirst){
2937 h[nJets] = (TH2F*)h[ij]->Clone(Form(cMask,nJets));
2938 }
2939 else{
2940 h[nJets]->Add(h[ij]);
2941 }
2942 }
2943 else{
2944 Printf("%s not found",Form(cMask,ij));
2945 }
2946 }
2947
2948}
2949
2950void ScaleNevent(TH1* h,TFile *fIn,Int_t iCond,Int_t iMC,Int_t iCen){
2951 // fetch the trigger histos
2952
2953 TDirectory *dir = (TDirectory*)fIn->Get("PWG4_services");
2954 TList *list1 = (TList*)dir->Get("pwg4serv");
2955
2956 TH2F* hTriggerCount =(TH2F*)list1->FindObject("fh2TriggerCount");
2957 Float_t nevents = hTriggerCount->GetBinContent(hTriggerCount->FindBin(0+iCen,iCond)); // we have to scale with number of triggers?
2958
2959 Float_t xsection = 1;
2960 if(iMC==7000){
2961 // 7 TeV sigma ND = 48.85 sigma NEL = 71.39 (Pythia 8... confirm with 6
2962 xsection = 71.39;
2963 Printf("MC Scaling setting nevents=%f to 1, xsection = %E",nevents,xsection);
2964 nevents = 1;
2965 }
2966
2967 Printf("Scaling %s with number of event %f",h->GetName(),nevents);
2968 if(nevents){
2969 Float_t scalef = 1./nevents/xsection;
2970 Printf("scale factor %E",scalef);
2971 h->Scale(scalef);
2972 }
2973
2974}
2975
c08c3ad2 2976Double_t Gamma0(Double_t *x,Double_t *par){
2977
2978 // Fixed multiplicity M
2979
2980
2981 Double_t p = par[0];
2982 Double_t b = par[1];
2983 Double_t A = par[2];
2984 Double_t xval = x[0] + p/b;
2985 // xval += par[3];
2986 if(xval<0)return 0;
2987
2988 // take log to avoid zeros and NANs
2989 Double_t f1 = TMath::Log(A*b);
8dc333da 2990 Double_t f11 = 1;// ROOT::Math::lgamma(p);
c08c3ad2 2991 Double_t f2 = TMath::Log(b * xval)*(p-1);
2992 Double_t f3 = -1.*b*xval;
2993 Double_t f = f1-f11+f2+f3;
2994 f = TMath::Exp(f);
2995 return f;
2996}
2997
2998Double_t Gamma(Double_t *x,Double_t *par){
2999
3000 // Fixed multiplicity M
3001 Double_t M = par[0];
3002 // Paramter p = <p_T>^2/sigma_p_T^2 adn b = <p_T>/sigma_p_T^2
3003 Double_t b = par[1]/par[2]/par[2];
3004 Double_t p = par[1] * b;
3005
8dc333da 3006 Double_t f = M * b / 1; // ROOT::Math::tgamma(M*p) * TMath::Power(M * b * x[0],M*(p-1)) * TMath::Exp(M*b*x[0]);
3007
3008 return f;
3009}
3010
3011
3012TF1* FitLHSgaus(TH1D *hDeltaPt, double minPt, double maxPt, int minIterations, int maxIterations, double maxDeltaMean, int mode, double minBoundSigma, double maxBoundSigma)
3013{
3014
3015 Double_t minPtBound = minPt;
3016 Double_t maxPtBound = maxPt;
3017
3018 Int_t nIter = 0;
3019 Double_t deltaMean = 999.;
3020 Double_t oldMean = 0.;
3021
3022 TF1 *f1 = new TF1("f1","gaus",minPtBound,maxPtBound);
3023 Double_t sigma = 0.;
3024 Double_t sigma_error = 0.;
3025 Double_t mean = 0.;
3026 Double_t mean_error = 0.;
3027
3028 while(nIter<minIterations || (nIter<maxIterations && deltaMean>maxDeltaMean)){
3029
3030 if(nIter>0){ // for initial fit use minPt and maxPt
3031 if(mode==0){
3032 maxPtBound = mean+5.;
3033 }
3034 if(mode==1){
3035 minPtBound = mean-(minBoundSigma*sigma);
3036 maxPtBound = mean+(maxBoundSigma*sigma);
3037 }
3038 }
3039
3040 f1->SetRange(minPtBound, maxPtBound);
3041 hDeltaPt->Fit(f1,"R0");
3042 Printf("fit range: %2.2f - %2.2f", minPtBound, maxPtBound);
3043 mean = f1->GetParameter(1);
3044 sigma = f1->GetParameter(2);
3045 mean_error = f1->GetParError(1);
3046 sigma_error = f1->GetParError(2);
3047
3048 deltaMean = TMath::Abs(mean-oldMean);
3049 oldMean = mean;
3050 Printf("FIT %d: #mu = %2.2f +- %2.2f (diff %2.2f), #sigma = %2.2f +- %2.2f, #chi_{2}/NDF = %2.2f/%d", nIter, mean, mean_error, deltaMean, sigma, sigma_error, f1->GetChisquare(), f1->GetNDF());
3051 nIter++;
3052 }
3053
3054 f1->SetRange(f1->GetX(1E-5,-70.,0.)-0.75,f1->GetX(1E-5,0.,70.)+0.75);
c08c3ad2 3055
8dc333da 3056 return f1;
c08c3ad2 3057}
3058
8dc333da 3059
3060
814a52a0 3061TF1* FitLHSgaus(TH1D *hDeltaPt, double &sigma, double &sigma_error, double &mean)
3062{
407d06f5 3063
3064 TF1 *f1 = new TF1("LHSgaus","gaus");
3065 f1->SetParameters(1,0,11);
3066 hDeltaPt->Fit(f1,"R0","",-50.,0.);
3067 // f1 = hDeltaPt->GetFunction("gaus");
814a52a0 3068 sigma = f1->GetParameter(2);
3069 mean = f1->GetParameter(1);
407d06f5 3070 hDeltaPt->Fit(f1,"R0","",mean-3.*sigma,mean+0.3*sigma);
814a52a0 3071 mean = f1->GetParameter(1);
3072 sigma = f1->GetParameter(2);
3073 sigma_error = f1->GetParError(2);
3074
3075 return f1;
3076}
8dc333da 3077
3078
3079
3080
3081TH2F* GetTH2PlotB(const char *cPath,Int_t embType, Int_t classType, Int_t cl, Int_t rp){
3082
3083 // emb type 0: single tracks, 1: emb jets
3084 // class type 0: centrality, 1: multiplicity (nb. of input tracks)
3085 // cl: centrality or multplicity class, -1 all
3086 // rp -1: all, 0: in plane, 1: in between, 2: out of plane
3087
3088 TString sClType;
3089 if(classType) sClType = "mult";
3090 else sClType = "cent";
3091 TString sEmbType;
3092 if(embType==0) sEmbType = "singles";
3093 if(embType==1) sEmbType = "jets";
3094
3095 TString path(cPath);
3096 TString file = Form("%s/jetresponse_%s_%s.root",path.Data(),sClType.Data(),sEmbType.Data());
3097 TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject(file.Data());
3098 if(!f)f = TFile::Open(file);
3099
3100 if(!f){
3101 Printf("GetTH2PlotB:%d Could not open file %s",__LINE__,file.Data());
3102 return 0x0;
3103 }
3104
3105 TString dfName = Form("JetsDeltaPt/rp %d/%s %d",rp,sClType.Data(),cl);
3106 TDirectory* df = (TDirectory*)f->Get(dfName.Data());
3107 if(!df){
3108 Printf("GetTH2PlotB:%d Could not open directory %s",__LINE__,dfName.Data());
3109 return 0x0;
3110 }
3111
3112 TString h2Name = "delta pT vs Probe pT";
3113 TH2F* h2 = (TH2F*)(df->Get(h2Name.Data()))->Clone(Form("h2B_%d",cl));
3114 if(!h2){
3115 Printf("GetTH2PlotB:%d Could not get histogram %s",__LINE__,h2Name.Data());
3116 return 0x0;
3117 }
3118
3119 return h2;
3120}
3121
8f206bf0 3122TObject* GetObjectFromFile(const char *cName,const char *cFile,const char* cDir,const char *cRep2,const char* cRep1){
8dc333da 3123 TDirectory *opwd = gDirectory;
3124 TFile *fIn = (TFile*)gROOT->GetListOfFiles()->FindObject(cFile);
3125 if(!fIn)fIn = TFile::Open(cFile);
3126 opwd->cd();
3127
3128 if(!fIn){
3129 Printf("File %s not found",cFile);
3130 return 0;
3131 }
3132 TDirectory *dir = (TDirectory*)fIn->Get(cDir);
3133 if(!dir){
3134 Printf("dir %s not found",cDir);
3135 return 0;
3136 }
3137 TString sList(cDir);
8f206bf0 3138 sList.ReplaceAll(cRep1,cRep2);
8dc333da 3139 TList *list = (TList*)dir->Get(sList.Data());
3140 if(!list){
3141 Printf("list %s not found",sList.Data());
3142 return 0;
3143 }
3144
3145 TObject *obj = list->FindObject(cName);
3146 if(!obj){
3147 Printf("object %s not found",cName);
3148 return 0;
3149 }
3150 return obj;
3151 }
3152
3153Double_t GetPoissonFluctuation(TH1 *h1,Double_t areaIn,Double_t areaJet){
3154 if(!h1)return 0;
3155 Printf(">>>> %s ",h1->GetName());
3156 Double_t meanPt = h1->GetMean();
3157 Double_t rmsPt = h1->GetRMS();
3158 Double_t mult = h1->Integral("width");
3159
3160 Double_t multJet = mult/areaIn*areaJet;
3161 Double_t sigma = TMath::Sqrt(multJet) * TMath::Sqrt(meanPt*meanPt+rmsPt*rmsPt);
3162 Printf("MeanPt %6.3f RMS %6.3f Tracks: %d",meanPt,rmsPt,(Int_t)mult);
3163 Printf("Tracks/jet: %d SigmaJet: %6.3f",(Int_t)multJet,sigma);
3164 return sigma;
3165}
3166
8f206bf0 3167void SetHistoAttributes(TH1* h1,Int_t iMarker,Int_t iColor,Float_t fMarkerSize){
8dc333da 3168 if(!h1)return;
8f206bf0 3169
3170
3171 h1->SetLineWidth(gStyle->GetLineWidth());
8dc333da 3172 h1->SetLineColor(iColor);
8f206bf0 3173
3174 h1->SetMarkerColor(iColor);
8dc333da 3175 h1->SetMarkerStyle(iMarker);
8f206bf0 3176 if(fMarkerSize>0){
3177 h1->SetMarkerSize(fMarkerSize);
3178 if(iMarker==29||iMarker==33) h1->SetMarkerSize(fMarkerSize * 1.2); // st
3179 }
8dc333da 3180}
3181
3182void SetGraphAttributes(TGraph* gr,Int_t iMarker,Int_t iColor){
3183 if(!gr)return;
8f206bf0 3184
8dc333da 3185 gr->SetLineColor(iColor);
8f206bf0 3186 gr->SetLineWidth(gStyle->GetLineWidth());
3187
3188 gr->SetMarkerColor(iColor);
8dc333da 3189 gr->SetMarkerStyle(iMarker);
8f206bf0 3190 gr->SetMarkerSize(1.7);
3191 if(iMarker==29||iMarker==33||iMarker==31) gr->SetMarkerSize(2.2); // st
8dc333da 3192}
3193
3194
3195void CloseFiles(){
3196 TSeqCollection *coll = gROOT->GetListOfFiles();
3197 for(int i = 0;i<coll->GetEntries();i++){
3198 TFile *f = (TFile*)coll->At(i);
3199 Printf("Closing %d %s",i,f->GetName());
3200 f->Close();
3201 }
3202}
8f206bf0 3203
3204
3205void SetHistAxisStyle(TH1* h1,Float_t titleOffsetX,Float_t titleSizeX,Float_t labelSizeX,
3206 Float_t titleOffsetY,Float_t titleSizeY,Float_t labelSizeY){
3207
3208
3209 if(titleOffsetX>0)h1->SetTitleOffset(titleOffsetX,"Y");
3210 if(titleOffsetY>0)h1->SetTitleOffset(titleOffsetY,"X");
3211
3212 if(titleSizeX>0)h1->SetTitleSize(h1->GetTitleSize("X")*titleSizeX,"X");
3213 if(titleSizeY>0)h1->SetTitleSize(h1->GetTitleSize("Y")*titleSizeX,"Y");
3214
3215 if(labelSizeX>0)h1->SetLabelSize(h1->GetLabelSize("X")*labelSizeX,"X");
3216 if(labelSizeY>0)h1->SetLabelSize(h1->GetLabelSize("Y")*labelSizeX,"Y");
3217
3218}
3219
3220void DrawDate(){
5f5e3e77 3221 static TDatime *time = new TDatime();
8f206bf0 3222 static TLatex *labelDate= new TLatex();
3223 labelDate->SetNDC(kTRUE);
3224 labelDate->SetTextAlign(11);
3225 labelDate->SetTextSize(0.04);
3226 labelDate->SetTextFont(52);
3227 labelDate->DrawLatex(0.01,0.01,Form("%d-%d",time->GetDate(),time->GetTime()));
3228}
3229
3230
5f5e3e77 3231void DrawALICE(TCanvas *c,Int_t iType,Float_t xCenter,Float_t yCenter,TString cExtra2,bool bExtraBox){
3232
3233 if(iType>10)DrawDate();
3234 Float_t aspectRatio = c->GetXsizeReal()/c->GetYsizeReal();;
3235
3236 /*
3237 if(bExtraBox){
3238 const float yHWB = 3.* 0.04;
3239 const float xHWB = 3.* 0.04/aspectRatio;
3240 TPad *p1 = new TPad("extraPad","extra",xCenter-xHWB, yCenter-yHWB,
3241 xCenter+xHWB, yCenter+yHWB);
3242 p1->SetFillColor(kBlue);
3243 p1->SetLineColor(kBlue);
3244 p1->Draw();
3245 }
3246 */
8f206bf0 3247 static TDatime *time = new TDatime();
3248 // iType = 1 Work in p
3249 //
3250 TString cType = "";
3251 TString cExtra1 = "";
3252 if((iType%10)==1){
3253 cType = "ALICE work in progress";
3254 }
3255 else if((iType%10)==2){
5f5e3e77 3256 cType = "ALICE performance";
8f206bf0 3257 cExtra1 = "Uncorrected";
3258 if(cExtra2.Length()==0){
3259 cExtra2 = Form("%04d/%02d/%02d",time->GetYear(),time->GetMonth(),time->GetDay());
3260 }
3261 }
5f5e3e77 3262 else if((iType%10)==3){
3263 cType = "ALICE preliminary";
8f206bf0 3264 }
3265 else{
3266 Printf("No ALICE type %d",iType);
3267 return;
3268 }
3269
3270
3271
3272 if(iType>10){
5f5e3e77 3273 cType = Form("%s requested",cType.Data());
8f206bf0 3274 }
3275
3276 if(cExtra2.Length()!=0&&cExtra1.Length()==0){
3277 cExtra1 = cExtra2;
3278 cExtra2 = "";
3279 }
3280
3281
5f5e3e77 3282 Float_t aspectRatio = c->GetXsizeReal()/c->GetYsizeReal();;
3283
3284 const float yHW = 1.5* 0.04;
3285 const float xHW = 1.5* 0.04/aspectRatio;
8f206bf0 3286
3287 TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",
5f5e3e77 3288 xCenter-xHW,yCenter-yHW,
3289