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