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