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