PWGLF/totEt/macros/emEt/PlotProtons.C

   1 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name);
   2 void SetStyles(TH1 *histo,int marker, int color){
   3   histo->Sumw2();
   4   histo->SetMarkerStyle(marker);
   5   histo->SetMarkerColor(color);
   6   histo->SetLineColor(color);
   7   //histo->GetXaxis()->SetTitle(xtitle);
   8   //histo->GetYaxis()->SetTitle(ytitle);
   9 }
  10 void Rescale(TH1 *histo,Int_t rescale){
  11   histo->Rebin(rescale);
  12   histo->Scale(1.0/((Float_t)rescale));
  13 }
  14 void SetStyles(TH1 *histo,int marker, int color,char *name){
  15   SetStyles(histo,marker,color);
  16   histo->SetName(name);
  17 }
  18 void PrintInfo(TH1 *histo){
  19   cout<<histo->GetName()<<":"<<endl;
  20   cout<<"x range "<<histo->GetBinLowEdge(1)<<" - "<<histo->GetBinLowEdge(histo->GetNbinsX()+1)<<endl;
  21   cout<<"y range "<<histo->GetMinimum()<<" - "<<histo->GetMaximum()<<endl;
  22 }
  23 Int_t colors[] = {0,TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack,
  24                     TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack,
  25                     TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack,
  26                     TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack};
  27 Int_t markers[] = {20,21,22,23,33, 24,25,26,32,27, 20,21,22,23,33, 24,25,26,32,27};
  28
  29
  30 Float_t nNeutronsSim[] = {0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0};
  31 Float_t nNeutronsData[] = {0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0};
  32 Float_t nNeutronsDataErr[] = {0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0};
  33 Float_t nNeutronsShortSim[] = {0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0};
  34 Float_t nNeutronsShortData[] = {0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0};
  35 Float_t nNeutronsShortDataErr[] = {0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0};
  36 Float_t nNeutronsSimCl[] = {0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0};
  37 Float_t nNeutronsDataCl[] = {0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0};
  38 Float_t nNeutronsShortSimCl[] = {0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0};
  39 Float_t nNeutronsShortDataCl[] = {0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0};
  40
  41 void WriteLatex();
  42 Float_t neutronCorrEmcal[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  43 Float_t neutronCorrPhos[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  44 Float_t neutronErrorEmcal[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  45 Float_t neutronErrorPhos[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  46 TH1D* Divide(TH1D* numerator, TH1D* denominator,Char_t* name){
  47   //TH1D *output = numerator->Clone(name);
  48   TH1D *output = new TH1D(name,name,numerator->GetNbinsX(),numerator->GetBinLowEdge(1),numerator->GetBinLowEdge(numerator->GetNbinsX()+1));
  49   output->GetXaxis()->SetTitle(numerator->GetXaxis()->GetTitle());
  50   output->GetYaxis()->SetTitle(numerator->GetYaxis()->GetTitle());
  51   Int_t nbinsNum = numerator->GetNbinsX();
  52   Int_t nbinsDen = denominator->GetNbinsX();
  53   for(Int_t i=1;i<=output->GetNbinsX();i++){
  54     if(nbinsNum!=nbinsDen){
  55       if(denominator->GetBinContent(i)>0){
  56         output->SetBinContent(i,numerator->GetBinContent(i)/denominator->GetBinContent(i));
  57       }
  58     }
  59     else{
  60       output->SetBinContent(i,numerator->GetBinContent(i));
  61     }
  62     //output->SetBinError((Int_t)i,0,0,(Double_t)0.0);
  63     //denominator->SetBinError((Int_t)i,0,0,(Double_t)0.0);
  64     //cout<<" "<<output->GetBinError(i)<<" "<<denominator->GetBinError(i)<<endl;
  65   }
  66   if(nbinsNum==nbinsDen){
  67     output->Divide(denominator);
  68   }
  69   //cout<<"Fixing "<<output->GetName()<<endl;
  70   Double_t newerror = 1e-5;
  71   for(Int_t i=1;i<=output->GetNbinsX();i++){
  72     //cout<<"Bin "<<i<<" setting bin error "<<output->GetBinError(i)<<" to 0.  New bin error: ";
  73     output->SetBinError((Int_t)i,(Double_t)newerror);
  74     //output->SetBinError((Int_t)i,0,0,(Double_t)0.0);
  75     //cout<<" "<<output->GetBinError(i)<<endl;
  76   }
  77   return output;
  78 }
  79
  80 void PlotProtons(TString filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root", TString filenameData = "rootFiles/LHC10hPass2/Et.ESD.realPbPb.EMCal.LHC10hPass2.Run139465.root", Bool_t effCorr = kTRUE){
  81   gStyle->SetOptTitle(0);
  82   gStyle->SetOptStat(0);
  83   gStyle->SetOptFit(0);
  84   Bool_t isPhos = kTRUE;
  85   TString detector = "Phos";
  86   if(filename.Contains("EMC")){
  87     detector = "Emcal";
  88     isPhos = kFALSE;
  89   }
  90
  91   TString tag = "";
  92   if(!effCorr) tag = "NoEffCorr";
  93
  94   ofstream myfile;
  95   TString textfilename = "Neutrons"+detector+tag+".dat";
  96   myfile.open (textfilename.Data());
  97   ofstream myfile2;
  98   TString textfilename2 = "Neutrons"+detector+tag+"Short.dat";
  99   myfile2.open (textfilename2.Data());
 100
 101   //Reading histograms
 102   TFile *f = TFile::Open(filename, "READ");
 103   TList *l = dynamic_cast<TList*>(f->Get("out1"));
 104   TH2F *fHistPIDProtonsTrackMatchedDepositedVsNch;
 105   TH2F *fHistPIDAntiProtonsTrackMatchedDepositedVsNch;
 106   TH2F *fHistPiKPTrackMatchedDepositedVsNch;
 107   TH2F *fHistNeutronsDepositedVsNch;
 108   TH2F *fHistAntiNeutronsDepositedVsNch;
 109   TH2F *fHistProtonsDepositedVsNch;
 110   TH2F *fHistAntiProtonsDepositedVsNch;
 111   TH2F *fHistPiKPDepositedVsNch;
 112   TH2F *fHistProtonsNotTrackMatchedDepositedVsNch;
 113   TH2F *fHistAntiProtonsNotTrackMatchedDepositedVsNch;
 114   TH2F *fHistPIDProtonsTrackMatchedDepositedVsNcl;
 115   TH2F *fHistPIDAntiProtonsTrackMatchedDepositedVsNcl;
 116   TH2F *fHistPiKPTrackMatchedDepositedVsNcl;
 117   TH2F *fHistNeutronsDepositedVsNcl;
 118   TH2F *fHistAntiNeutronsDepositedVsNcl;
 119   TH2F *fHistProtonsDepositedVsNcl;
 120   TH2F *fHistAntiProtonsDepositedVsNcl;
 121   TH2F *fHistPiKPDepositedVsNcl;
 122   TH2F *fHistProtonsNotTrackMatchedDepositedVsNcl;
 123   TH2F *fHistAntiProtonsNotTrackMatchedDepositedVsNcl;
 124
 125   if(effCorr){
 126     fHistPIDProtonsTrackMatchedDepositedVsNch = (TH2F *)l->FindObject("fHistPIDProtonsTrackMatchedDepositedVsNch");
 127     fHistPIDAntiProtonsTrackMatchedDepositedVsNch = (TH2F *)l->FindObject("fHistPIDAntiProtonsTrackMatchedDepositedVsNch");
 128     fHistPiKPTrackMatchedDepositedVsNch = (TH2F *)l->FindObject("fHistPiKPTrackMatchedDepositedVsNch");
 129     fHistNeutronsDepositedVsNch = (TH2F *)l->FindObject("fHistNeutronsDepositedVsNch");
 130     fHistAntiNeutronsDepositedVsNch = (TH2F *)l->FindObject("fHistAntiNeutronsDepositedVsNch");
 131     fHistProtonsDepositedVsNch = (TH2F *)l->FindObject("fHistProtonsDepositedVsNch");
 132     fHistAntiProtonsDepositedVsNch = (TH2F *)l->FindObject("fHistAntiProtonsDepositedVsNch");
 133     fHistPiKPDepositedVsNch = (TH2F *)l->FindObject("fHistPiKPDepositedVsNch");
 134     fHistProtonsNotTrackMatchedDepositedVsNch = (TH2F *)l->FindObject("fHistProtonsNotTrackMatchedDepositedVsNch");
 135     fHistAntiProtonsNotTrackMatchedDepositedVsNch = (TH2F *)l->FindObject("fHistAntiProtonsNotTrackMatchedDepositedVsNch");
 136
 137     fHistPIDProtonsTrackMatchedDepositedVsNcl = (TH2F *)l->FindObject("fHistPIDProtonsTrackMatchedDepositedVsNcl");
 138     fHistPIDAntiProtonsTrackMatchedDepositedVsNcl = (TH2F *)l->FindObject("fHistPIDAntiProtonsTrackMatchedDepositedVsNcl");
 139     fHistPiKPTrackMatchedDepositedVsNcl = (TH2F *)l->FindObject("fHistPiKPTrackMatchedDepositedVsNcl");
 140     fHistNeutronsDepositedVsNcl = (TH2F *)l->FindObject("fHistNeutronsDepositedVsNcl");
 141     fHistAntiNeutronsDepositedVsNcl = (TH2F *)l->FindObject("fHistAntiNeutronsDepositedVsNcl");
 142     fHistProtonsDepositedVsNcl = (TH2F *)l->FindObject("fHistProtonsDepositedVsNcl");
 143     fHistAntiProtonsDepositedVsNcl = (TH2F *)l->FindObject("fHistAntiProtonsDepositedVsNcl");
 144     fHistPiKPDepositedVsNcl = (TH2F *)l->FindObject("fHistPiKPDepositedVsNcl");
 145     fHistProtonsNotTrackMatchedDepositedVsNcl = (TH2F *)l->FindObject("fHistProtonsNotTrackMatchedDepositedVsNcl");
 146     fHistAntiProtonsNotTrackMatchedDepositedVsNcl = (TH2F *)l->FindObject("fHistAntiProtonsNotTrackMatchedDepositedVsNcl");
 147   }
 148   else{
 149     fHistPIDProtonsTrackMatchedDepositedVsNch = (TH2F *)l->FindObject("fHistPIDProtonsTrackMatchedDepositedVsNchNoEff");
 150     fHistPIDAntiProtonsTrackMatchedDepositedVsNch = (TH2F *)l->FindObject("fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff");
 151     fHistPiKPTrackMatchedDepositedVsNch = (TH2F *)l->FindObject("fHistPiKPTrackMatchedDepositedVsNchNoEff");
 152     fHistNeutronsDepositedVsNch = (TH2F *)l->FindObject("fHistNeutronsDepositedVsNchNoEffCorr");
 153     fHistAntiNeutronsDepositedVsNch = (TH2F *)l->FindObject("fHistAntiNeutronsDepositedVsNchNoEffCorr");
 154     fHistProtonsDepositedVsNch = (TH2F *)l->FindObject("fHistProtonsDepositedVsNchNoEffCorr");
 155     fHistAntiProtonsDepositedVsNch = (TH2F *)l->FindObject("fHistAntiProtonsDepositedVsNchNoEffCorr");
 156     fHistPiKPDepositedVsNch = (TH2F *)l->FindObject("fHistPiKPDepositedVsNchNoEffCorr");
 157     fHistProtonsNotTrackMatchedDepositedVsNch = (TH2F *)l->FindObject("fHistProtonsNotTrackMatchedDepositedVsNchNoEffCorr");
 158     fHistAntiProtonsNotTrackMatchedDepositedVsNch = (TH2F *)l->FindObject("fHistAntiProtonsNotTrackMatchedDepositedVsNchNoEffCorr");
 159
 160
 161
 162     fHistPIDProtonsTrackMatchedDepositedVsNcl = (TH2F *)l->FindObject("fHistPIDProtonsTrackMatchedDepositedVsNclNoEff");
 163     fHistPIDAntiProtonsTrackMatchedDepositedVsNcl = (TH2F *)l->FindObject("fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff");
 164     fHistPiKPTrackMatchedDepositedVsNcl = (TH2F *)l->FindObject("fHistPiKPTrackMatchedDepositedVsNclNoEff");
 165     fHistNeutronsDepositedVsNcl = (TH2F *)l->FindObject("fHistNeutronsDepositedVsNclNoEffCorr");
 166     fHistAntiNeutronsDepositedVsNcl = (TH2F *)l->FindObject("fHistAntiNeutronsDepositedVsNclNoEffCorr");
 167     fHistProtonsDepositedVsNcl = (TH2F *)l->FindObject("fHistProtonsDepositedVsNclNoEffCorr");
 168     fHistAntiProtonsDepositedVsNcl = (TH2F *)l->FindObject("fHistAntiProtonsDepositedVsNclNoEffCorr");
 169     fHistPiKPDepositedVsNcl = (TH2F *)l->FindObject("fHistPiKPDepositedVsNclNoEffCorr");
 170     fHistProtonsNotTrackMatchedDepositedVsNcl = (TH2F *)l->FindObject("fHistProtonsNotTrackMatchedDepositedVsNclNoEffCorr");
 171     fHistAntiProtonsNotTrackMatchedDepositedVsNcl = (TH2F *)l->FindObject("fHistAntiProtonsNotTrackMatchedDepositedVsNclNoEffCorr");
 172   }
 173
 174   TH3F *fHistCentVsNchVsNcl = l->FindObject("fHistCentVsNchVsNclReco");
 175   fHistCentVsNchVsNcl->GetXaxis()->SetTitle("cent");
 176   fHistCentVsNchVsNcl->GetYaxis()->SetTitle("N_{Ch}");
 177   fHistCentVsNchVsNcl->GetZaxis()->SetTitle("N_{Cl}");
 178   //fHistCentVsNchVsNcl->SetName("fHistCentVsNchVsNclSim");
 179
 180   fHistPIDProtonsTrackMatchedDepositedVsNch->SetName(Form("%sSim","fHistPIDProtonsTrackMatchedDepositedVsNch"));
 181   fHistPIDAntiProtonsTrackMatchedDepositedVsNch->SetName(Form("%sSim","fHistPIDAntiProtonsTrackMatchedDepositedVsNch"));
 182   fHistPiKPTrackMatchedDepositedVsNch->SetName(Form("%sSim","fHistPiKPTrackMatchedDepositedVsNch"));
 183
 184   TFile *fData = TFile::Open(filenameData, "READ");
 185   TList *lData = dynamic_cast<TList*>(fData->Get("out1"));
 186   TH2F *fHistPIDProtonsTrackMatchedDepositedVsNchData;
 187   TH2F *fHistPIDAntiProtonsTrackMatchedDepositedVsNchData;
 188   TH2F *fHistPiKPTrackMatchedDepositedVsNchData;
 189   TH2F *fHistPIDProtonsTrackMatchedDepositedVsNclData;
 190   TH2F *fHistPIDAntiProtonsTrackMatchedDepositedVsNclData;
 191   TH2F *fHistPiKPTrackMatchedDepositedVsNclData;
 192   if(effCorr){
 193     fHistPIDProtonsTrackMatchedDepositedVsNchData =(TH2F *) lData->FindObject("fHistPIDProtonsTrackMatchedDepositedVsNch");
 194     fHistPIDAntiProtonsTrackMatchedDepositedVsNchData = (TH2F *)lData->FindObject("fHistPIDAntiProtonsTrackMatchedDepositedVsNch");
 195     fHistPiKPTrackMatchedDepositedVsNchData = (TH2F *)lData->FindObject("fHistPiKPTrackMatchedDepositedVsNch");
 196     fHistPIDProtonsTrackMatchedDepositedVsNclData =(TH2F *) lData->FindObject("fHistPIDProtonsTrackMatchedDepositedVsNcl");
 197     fHistPIDAntiProtonsTrackMatchedDepositedVsNclData = (TH2F *)lData->FindObject("fHistPIDAntiProtonsTrackMatchedDepositedVsNcl");
 198     fHistPiKPTrackMatchedDepositedVsNclData = (TH2F *)lData->FindObject("fHistPiKPTrackMatchedDepositedVsNcl");
 199   }
 200   else{
 201     fHistPIDProtonsTrackMatchedDepositedVsNchData =(TH2F *) lData->FindObject("fHistPIDProtonsTrackMatchedDepositedVsNchNoEff");
 202     fHistPIDAntiProtonsTrackMatchedDepositedVsNchData = (TH2F *)lData->FindObject("fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff");
 203     fHistPiKPTrackMatchedDepositedVsNchData = (TH2F *)lData->FindObject("fHistPiKPTrackMatchedDepositedVsNchNoEff");
 204     fHistPIDProtonsTrackMatchedDepositedVsNclData =(TH2F *) lData->FindObject("fHistPIDProtonsTrackMatchedDepositedVsNclNoEff");
 205     fHistPIDAntiProtonsTrackMatchedDepositedVsNclData = (TH2F *)lData->FindObject("fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff");
 206     fHistPiKPTrackMatchedDepositedVsNclData = (TH2F *)lData->FindObject("fHistPiKPTrackMatchedDepositedVsNclNoEff");
 207   }
 208   TH3F *fHistCentVsNchVsNclData = lData->FindObject("fHistCentVsNchVsNclReco");
 209   fHistPIDProtonsTrackMatchedDepositedVsNchData->SetName(Form("%sData","fHistPIDProtonsTrackMatchedDepositedVsNch"));
 210   fHistPIDAntiProtonsTrackMatchedDepositedVsNchData->SetName(Form("%sData","fHistPIDAntiProtonsTrackMatchedDepositedVsNch"));
 211   fHistPiKPTrackMatchedDepositedVsNchData->SetName(Form("%sData","fHistPiKPTrackMatchedDepositedVsNch"));
 212   fHistPIDProtonsTrackMatchedDepositedVsNclData->SetName(Form("%sData","fHistPIDProtonsTrackMatchedDepositedVsNcl"));
 213   fHistPIDAntiProtonsTrackMatchedDepositedVsNclData->SetName(Form("%sData","fHistPIDAntiProtonsTrackMatchedDepositedVsNcl"));
 214   //fHistPiKPTrackMatchedDepositedVsNclData->SetName(Form("%sData","fHistPiKPTrackMatchedDepositedVsNcl"));
 215   fHistCentVsNchVsNclData->SetName("fHistCentVsNchVsNclRecoData");
 216   //TH3F *fHistCentVsNchVsNclData = lData->FindObject("fHistCentVsNchVsNcl");
 217
 218   TCanvas *c1 = new TCanvas("c1","Simulation",600,400);
 219   c1->SetTopMargin(0.02);
 220   c1->SetRightMargin(0.149329);
 221   c1->SetBorderSize(0);
 222   c1->SetFillColor(0);
 223   c1->SetFillColor(0);
 224   c1->SetBorderMode(0);
 225   c1->SetFrameFillColor(0);
 226   c1->SetFrameBorderMode(0);
 227   //c1->SetLogz();
 228   //fHistPIDProtonsTrackMatchedDepositedVsNch->Draw("colz");
 229   TH1D *fHistPiKPTrackMatchedDepositedVsNchProf = fHistPiKPTrackMatchedDepositedVsNch->ProfileY();
 230   fHistPiKPTrackMatchedDepositedVsNchProf->Scale(1.0/30.0);
 231   fHistPiKPTrackMatchedDepositedVsNchProf->GetXaxis()->SetTitle("N_{ch}");
 232   fHistPiKPTrackMatchedDepositedVsNchProf->GetYaxis()->SetTitle("<E_{T}>");
 233   TH1D *fHistPiKPTrackMatchedDepositedVsNchProfCopy = fHistPiKPTrackMatchedDepositedVsNchProf->Clone("fHistPiKPTrackMatchedDepositedVsNchProfCopy");
 234   fHistPiKPTrackMatchedDepositedVsNchProfCopy->Draw();
 235   TH1D *fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf = fHistPIDAntiProtonsTrackMatchedDepositedVsNch->ProfileY("fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf");
 236   TH1D *fHistPIDAntiProtonsTrackMatchedDepositedVsNclProf = fHistPIDAntiProtonsTrackMatchedDepositedVsNcl->ProfileY("fHistPIDAntiProtonsTrackMatchedDepositedVsNclProf");
 237   fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf->Draw("same");
 238   TH1D *fHistPIDProtonsTrackMatchedDepositedVsNchProf = fHistPIDProtonsTrackMatchedDepositedVsNch->ProfileY("fHistPIDProtonsTrackMatchedDepositedVsNchProf");
 239   TH1D *fHistPIDProtonsTrackMatchedDepositedVsNclProf = fHistPIDProtonsTrackMatchedDepositedVsNcl->ProfileY("fHistPIDProtonsTrackMatchedDepositedVsNclProf");
 240   fHistPIDProtonsTrackMatchedDepositedVsNchProf->Draw("same");
 241   SetStyles(fHistPiKPTrackMatchedDepositedVsNchProf,29,1);
 242   SetStyles(fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf,20,TColor::kRed);
 243   SetStyles(fHistPIDProtonsTrackMatchedDepositedVsNchProf,24,TColor::kBlue);
 244   TLegend *leg = new TLegend(0.157718,0.709677,0.278523,0.935484);
 245   leg->SetFillStyle(0);
 246   leg->SetFillColor(0);
 247   leg->SetBorderSize(0);
 248   leg->SetTextSize(0.03);
 249   leg->SetTextSize(0.038682);
 250   leg->AddEntry(fHistPiKPTrackMatchedDepositedVsNchProf,"Track matched E_{T}/30");
 251   leg->AddEntry(fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf,"Identified #bar{p} E_{T}");
 252   leg->AddEntry(fHistPIDProtonsTrackMatchedDepositedVsNchProf,"Identified p E_{T}");
 253   leg->Draw();
 254   TString name1 = "/tmp/Sim"+detector+".png";
 255   c1->SaveAs(name1.Data());
 256   //cerr<<"176"<<endl;
 257
 258   TCanvas *c2 = new TCanvas("c2","Data",600,400);
 259   c2->SetTopMargin(0.02);
 260   c2->SetRightMargin(0.149329);
 261   c2->SetBorderSize(0);
 262   c2->SetFillColor(0);
 263   c2->SetFillColor(0);
 264   c2->SetBorderMode(0);
 265   c2->SetFrameFillColor(0);
 266   c2->SetFrameBorderMode(0);
 267   //c2->SetLogz();
 268   //fHistPIDProtonsTrackMatchedDepositedVsNchData->Draw("colz");
 269   TH1D *fHistPiKPTrackMatchedDepositedVsNchDataProf = fHistPiKPTrackMatchedDepositedVsNchData->ProfileY("fHistPiKPTrackMatchedDepositedVsNchDataProf");
 270   fHistPiKPTrackMatchedDepositedVsNchDataProf->Scale(1.0/30.0);
 271   fHistPiKPTrackMatchedDepositedVsNchDataProf->GetXaxis()->SetTitle("N_{ch}");
 272   fHistPiKPTrackMatchedDepositedVsNchDataProf->GetYaxis()->SetTitle("<E_{T}>");
 273   TH1D *fHistPiKPTrackMatchedDepositedVsNchDataProfCopy = fHistPiKPTrackMatchedDepositedVsNchDataProf->Clone("fHistPiKPTrackMatchedDepositedVsNchDataProfCopy");
 274   fHistPiKPTrackMatchedDepositedVsNchDataProfCopy->Draw();
 275   TH1D *fHistPIDAntiProtonsTrackMatchedDepositedVsNchDataProf = fHistPIDAntiProtonsTrackMatchedDepositedVsNchData->ProfileY("fHistPIDAntiProtonsTrackMatchedDepositedVsNchDataProf");
 276   TH1D *fHistPIDAntiProtonsTrackMatchedDepositedVsNclDataProf = fHistPIDAntiProtonsTrackMatchedDepositedVsNclData->ProfileY("fHistPIDAntiProtonsTrackMatchedDepositedVsNclDataProf");
 277   fHistPIDAntiProtonsTrackMatchedDepositedVsNchDataProf->Draw("same");
 278   TH1D *fHistPIDProtonsTrackMatchedDepositedVsNchDataProf = fHistPIDProtonsTrackMatchedDepositedVsNchData->ProfileY("fHistPIDProtonsTrackMatchedDepositedVsNchDataProf");
 279   TH1D *fHistPIDProtonsTrackMatchedDepositedVsNclDataProf = fHistPIDProtonsTrackMatchedDepositedVsNclData->ProfileY("fHistPIDProtonsTrackMatchedDepositedVsNclDataProf");
 280   fHistPIDProtonsTrackMatchedDepositedVsNchDataProf->Draw("same");
 281   SetStyles(fHistPiKPTrackMatchedDepositedVsNchDataProf,29,1);
 282   SetStyles(fHistPIDAntiProtonsTrackMatchedDepositedVsNchDataProf,20,TColor::kRed);
 283   SetStyles(fHistPIDProtonsTrackMatchedDepositedVsNchDataProf,24,TColor::kBlue);
 284   TLegend *legData = new TLegend(0.157718,0.709677,0.278523,0.935484);
 285   legData->SetFillStyle(0);
 286   legData->SetFillColor(0);
 287   legData->SetBorderSize(0);
 288   legData->SetTextSize(0.03);
 289   legData->SetTextSize(0.038682);
 290   legData->AddEntry(fHistPiKPTrackMatchedDepositedVsNchDataProf,"Track matched E_{T}/30");
 291   legData->AddEntry(fHistPIDAntiProtonsTrackMatchedDepositedVsNchDataProf,"Identified #bar{p} E_{T}");
 292   legData->AddEntry(fHistPIDProtonsTrackMatchedDepositedVsNchDataProf,"Identified p E_{T}");
 293   legData->Draw();
 294   //cerr<<"212"<<endl;
 295
 296   TString name2 = "/tmp/Data"+detector+".png";
 297   c2->SaveAs(name2.Data());
 298   fHistPiKPTrackMatchedDepositedVsNchDataProf->Scale(30.0);
 299   fHistPiKPTrackMatchedDepositedVsNchProf->Scale(30.0);
 300
 301   TCanvas *c3 = new TCanvas("c3","Ratios at low pT",600,400);
 302   c3->SetTopMargin(0.02);
 303   c3->SetRightMargin(0.149329);
 304   c3->SetBorderSize(0);
 305   c3->SetFillColor(0);
 306   c3->SetFillColor(0);
 307   c3->SetBorderMode(0);
 308   c3->SetFrameFillColor(0);
 309   c3->SetFrameBorderMode(0);
 310   TH1D *hAntiProtonSim = Divide(fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf,fHistPiKPTrackMatchedDepositedVsNchProf,"hAntiProtonSim");
 311   TH1D *hProtonSim = Divide(fHistPIDProtonsTrackMatchedDepositedVsNchProf,fHistPiKPTrackMatchedDepositedVsNchProf,"hProtonSim");
 312   TH1D *htemp =  fHistPIDProtonsTrackMatchedDepositedVsNchProf->Clone("temp");
 313   //cout<<htemp->GetEntries();
 314   htemp->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf);
 315   htemp->Scale(2);//because the profile averages when you add two profiles
 316   //cout<<" "<<htemp->GetEntries()<<endl;
 317   //cout<<"entries "<<fHistPiKPTrackMatchedDepositedVsNchProf->GetEntries()<<endl;
 318   TH1D *hAllProtonSim = Divide(htemp,fHistPiKPTrackMatchedDepositedVsNchProf,"hAllProtonSim");
 319   //delete htemp;
 320
 321
 322   //cerr<<"240"<<endl;
 323
 324   TH1D *hAntiProtonData = Divide(fHistPIDAntiProtonsTrackMatchedDepositedVsNchDataProf,fHistPiKPTrackMatchedDepositedVsNchDataProf,"hAntiProtonData");
 325   TH1D *hProtonData = Divide(fHistPIDProtonsTrackMatchedDepositedVsNchDataProf,fHistPiKPTrackMatchedDepositedVsNchDataProf,"hProtonData");
 326   TH1D *htemp2 =  fHistPIDProtonsTrackMatchedDepositedVsNchDataProf->Clone("temp2");
 327
 328   //cout<<htemp2->GetEntries();
 329   htemp2->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNchDataProf);
 330   htemp2->Scale(2);//because the profile averages when you add two profiles
 331   //cout<<" "<<htemp2->GetEntries()<<endl;
 332   //cout<<"entries "<<fHistPiKPTrackMatchedDepositedVsNchDataProf->GetEntries()<<endl;
 333   TH1D *hAllProtonData = Divide(htemp2,fHistPiKPTrackMatchedDepositedVsNchDataProf,"hAllProtonData");
 334   //delete htemp2;
 335   SetStyles(hAntiProtonSim,20,TColor::kRed);
 336   SetStyles(hProtonSim,24,TColor::kRed);
 337   SetStyles(hAllProtonSim,28,1);
 338   SetStyles(hAntiProtonData,22,TColor::kRed);
 339   SetStyles(hProtonData,26,TColor::kRed);
 340   SetStyles(hAllProtonData,33,1);
 341   hAllProtonSim->GetXaxis()->SetTitle("N_ch");
 342   hAllProtonSim->GetYaxis()->SetTitle("E_{T}^{p,#bar{p}}/E_{T}^{#pi,K,p TM}");
 343   TF1 *funcData = new TF1("funcData","[0]*exp(-x/[1])+[2]-[3]*x",0,2500);
 344   funcData->SetParameter(0,0.12);
 345   funcData->SetParameter(1,10);
 346   funcData->SetParLimits(1,1e-9,1e9);
 347   funcData->SetParameter(2,0.09);
 348   funcData->SetParameter(3,2e-5);
 349   funcData->SetLineColor(hAllProtonData->GetLineColor());
 350   hAllProtonData->Fit(funcData,"","",50,2700);
 351   TF1 *funcSim = new TF1("funcSim","[0]*exp(-x/[1])+[2]-[3]*x",0,2500);
 352   funcSim->SetParameter(0,0.06);
 353   funcSim->SetParLimits(0,0,0.12);
 354   funcSim->SetParameter(1,10);
 355   funcSim->SetParLimits(1,1e-3,1e2);
 356   funcSim->SetParameter(2,0.09);
 357   funcSim->SetParameter(3,2e-5);
 358   funcSim->SetLineColor(hAllProtonSim->GetLineColor());
 359   hAllProtonSim->Fit(funcSim,"","",0,2700);
 360   hAllProtonSim->Draw();
 361   hAllProtonData->Draw("same");
 362   hAntiProtonData->Draw("same");
 363   hProtonData->Draw("same");
 364   hAntiProtonSim->Draw("same");
 365   hProtonSim->Draw("same");
 366   TLegend *legRatio = new TLegend(0.157718,0.709677,0.278523,0.935484);
 367   legRatio->SetFillStyle(0);
 368   legRatio->SetFillColor(0);
 369   legRatio->SetBorderSize(0);
 370   legRatio->SetTextSize(0.03);
 371   legRatio->SetTextSize(0.038682);
 372   legRatio->AddEntry(hAllProtonSim,"p+#bar{p} Sim");
 373   legRatio->AddEntry(hAllProtonData,"p+#bar{p} Data");
 374   legRatio->AddEntry(hProtonSim,"p Sim");
 375   legRatio->AddEntry(hProtonData,"p Data");
 376   legRatio->AddEntry(hAntiProtonSim,"#bar{p} Sim");
 377   legRatio->AddEntry(hAntiProtonData,"#bar{p} Data");
 378   legRatio->Draw();
 379   //     c1->cd();
 380   //     SetStyles(htemp,21,1);
 381   //     htemp->Draw("same");
 382   //     c2->cd();
 383   //     SetStyles(htemp2,21,1);
 384   //     htemp2->Draw("same");
 385
 386   TString name3 = "/tmp/ProtonRatio"+detector+".eps";
 387   c3->SaveAs(name3.Data());
 388
 389   cerr<<"306"<<endl;
 390
 391
 392   TCanvas *c4 = new TCanvas("c4","Sim: ratios of particles to low pT protons",600,400);
 393   c4->SetTopMargin(0.02);
 394   c4->SetRightMargin(0.149329);
 395   c4->SetBorderSize(0);
 396   c4->SetFillColor(0);
 397   c4->SetFillColor(0);
 398   c4->SetBorderMode(0);
 399   c4->SetFrameFillColor(0);
 400   c4->SetFrameBorderMode(0);
 401   TH1D *fHistNeutronsDepositedVsNchProf = fHistNeutronsDepositedVsNch->ProfileY("fHistNeutronsDepositedVsNchProf");
 402   TH1D *fHistAntiNeutronsDepositedVsNchProf = fHistAntiNeutronsDepositedVsNch->ProfileY("fHistAntiNeutronsDepositedVsNchProf");
 403   TH1D *fHistProtonsDepositedVsNchProf = fHistProtonsDepositedVsNch->ProfileY("fHistProtonsDepositedVsNchProf");
 404   TH1D *fHistAntiProtonsDepositedVsNchProf = fHistAntiProtonsDepositedVsNch->ProfileY("fHistAntiProtonsDepositedVsNchProf");
 405   TH1D *hNeutronsOverLowPtProtons = Divide(fHistNeutronsDepositedVsNchProf,fHistPIDProtonsTrackMatchedDepositedVsNchProf,"hAllNeutronsOverLowPtProtonsName");
 406   TH1D *hAntiNeutronsOverLowPtAntiProtons = Divide(fHistAntiNeutronsDepositedVsNchProf,fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf,"hAllNeutronsOverLowPtAntiProtonsName");
 407   TH1D *hProtonsOverLowPtProtons = Divide(fHistProtonsDepositedVsNchProf,fHistPIDProtonsTrackMatchedDepositedVsNchProf,"hAllProtonsOverLowPtProtonsName");
 408   TH1D *hAntiProtonsOverLowPtAntiProtons = Divide(fHistAntiProtonsDepositedVsNchProf,fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf,"hAllProtonsOverLowPtAntiProtonsName");
 409
 410   TH1D *fHistNeutronsDepositedVsNclProf = fHistNeutronsDepositedVsNcl->ProfileY("fHistNeutronsDepositedVsNclProf");
 411   TH1D *fHistAntiNeutronsDepositedVsNclProf = fHistAntiNeutronsDepositedVsNcl->ProfileY("fHistAntiNeutronsDepositedVsNclProf");
 412   TH1D *fHistProtonsDepositedVsNclProf = fHistProtonsDepositedVsNcl->ProfileY("fHistProtonsDepositedVsNclProf");
 413   TH1D *fHistAntiProtonsDepositedVsNclProf = fHistAntiProtonsDepositedVsNcl->ProfileY("fHistAntiProtonsDepositedVsNclProf");
 414   cout<<"Dividing cluster histos"<<endl;
 415   TH1D *hNeutronsOverLowPtProtonsCl = Divide(fHistNeutronsDepositedVsNclProf,fHistPIDProtonsTrackMatchedDepositedVsNclProf,"hAllNeutronsOverLowPtProtonsClName");
 416   TH1D *hAntiNeutronsOverLowPtAntiProtonsCl = Divide(fHistAntiNeutronsDepositedVsNclProf,fHistPIDAntiProtonsTrackMatchedDepositedVsNclProf,"hAllNeutronsOverLowPtAntiProtonsClName");
 417   TH1D *hProtonsOverLowPtProtonsCl = Divide(fHistProtonsDepositedVsNclProf,fHistPIDProtonsTrackMatchedDepositedVsNclProf,"hAllProtonsOverLowPtProtonsClName");
 418   TH1D *hAntiProtonsOverLowPtAntiProtonsCl = Divide(fHistAntiProtonsDepositedVsNclProf,fHistPIDAntiProtonsTrackMatchedDepositedVsNclProf,"hAllProtonsOverLowPtAntiProtonsClName");
 419   cout<<"Done dividing cluster histos"<<endl;
 420
 421   TH1D *hTmpAllLowPtProtons = fHistPIDProtonsTrackMatchedDepositedVsNchProf->Clone("hTmpAllLowPtProtons");
 422   hTmpAllLowPtProtons->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf);
 423   TH1D *hTmpAllLowPtProtonsData = fHistPIDProtonsTrackMatchedDepositedVsNchDataProf->Clone("hTmpAllLowPtProtonsData");
 424   hTmpAllLowPtProtonsData->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNchDataProf);
 425   TH1D *hTmpAllProtons = fHistProtonsDepositedVsNchProf->Clone("hTmpAllProtons");
 426   hTmpAllProtons->Add(fHistAntiProtonsDepositedVsNchProf);
 427   TH1D *hTmpAllNeutrons = fHistNeutronsDepositedVsNchProf->Clone("hTmpAllNeutrons");
 428   hTmpAllNeutrons->Add(fHistAntiNeutronsDepositedVsNchProf);
 429
 430   TH1D *hTmpAllLowPtProtonsCl = fHistPIDProtonsTrackMatchedDepositedVsNclProf->Clone("hTmpAllLowPtProtonsCl");
 431   hTmpAllLowPtProtonsCl->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNclProf);
 432   TH1D *hTmpAllLowPtProtonsDataCl = fHistPIDProtonsTrackMatchedDepositedVsNclDataProf->Clone("hTmpAllLowPtProtonsDataCl");
 433   hTmpAllLowPtProtonsDataCl->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNclDataProf);
 434   TH1D *hTmpAllProtonsCl = fHistProtonsDepositedVsNclProf->Clone("hTmpAllProtonsCl");
 435   hTmpAllProtonsCl->Add(fHistAntiProtonsDepositedVsNclProf);
 436   TH1D *hTmpAllNeutronsCl = fHistNeutronsDepositedVsNclProf->Clone("hTmpAllNeutronsCl");
 437   hTmpAllNeutronsCl->Add(fHistAntiNeutronsDepositedVsNclProf);
 438
 439   TH1D *hAllNeutronsOverLowPtProtons = Divide(hTmpAllNeutrons,hTmpAllLowPtProtons,"hAllNeutronsOverLowPtProtons");
 440   TH1D *hAllNeutronsOverLowPtProtonsData = Divide(hTmpAllNeutrons,hTmpAllLowPtProtonsData,"hAllNeutronsOverLowPtProtonsData");
 441   TH1D *hAllProtonsOverLowPtProtons = Divide(hTmpAllProtons,hTmpAllLowPtProtons,"hAllProtonsOverLowPtProtons");
 442   TH1D *hAllProtonsOverLowPtProtonsData = Divide(hTmpAllProtons,hTmpAllLowPtProtonsData,"hAllProtonsOverLowPtProtonsData");
 443   SetStyles(hAntiNeutronsOverLowPtAntiProtons,20,TColor::kBlue);
 444   SetStyles(hNeutronsOverLowPtProtons,24,TColor::kBlue);
 445   SetStyles(hAntiProtonsOverLowPtAntiProtons,22,TColor::kRed);
 446   SetStyles(hProtonsOverLowPtProtons,26,TColor::kRed);
 447   SetStyles(hAllNeutronsOverLowPtProtonsData,29,TColor::kGreen+3);
 448   SetStyles(hAllNeutronsOverLowPtProtons,30,TColor::kGreen+3);
 449   SetStyles(hAllProtonsOverLowPtProtons,24,TColor::kRed);
 450   SetStyles(hAllProtonsOverLowPtProtonsData,20,TColor::kRed);
 451
 452   TH1D *hAllNeutronsOverLowPtProtonsCl = Divide(hTmpAllNeutronsCl,hTmpAllLowPtProtonsCl,"hAllNeutronsOverLowPtProtonsCl");
 453   TH1D *hAllNeutronsOverLowPtProtonsDataCl = Divide(hTmpAllNeutronsCl,hTmpAllLowPtProtonsDataCl,"hAllNeutronsOverLowPtProtonsDataCl");
 454   TH1D *hAllProtonsOverLowPtProtonsCl = Divide(hTmpAllProtonsCl,hTmpAllLowPtProtonsCl,"hAllProtonsOverLowPtProtonsCl");
 455   TH1D *hAllProtonsOverLowPtProtonsDataCl = Divide(hTmpAllProtonsCl,hTmpAllLowPtProtonsDataCl,"hAllProtonsOverLowPtProtonsDataCl");
 456   SetStyles(hAntiNeutronsOverLowPtAntiProtonsCl,20,TColor::kBlue);
 457   SetStyles(hNeutronsOverLowPtProtonsCl,24,TColor::kBlue);
 458   SetStyles(hAntiProtonsOverLowPtAntiProtonsCl,22,TColor::kRed);
 459   SetStyles(hProtonsOverLowPtProtonsCl,26,TColor::kRed);
 460   SetStyles(hAllNeutronsOverLowPtProtonsDataCl,29,TColor::kGreen+3);
 461   SetStyles(hAllNeutronsOverLowPtProtonsCl,30,TColor::kGreen+3);
 462   SetStyles(hAllProtonsOverLowPtProtonsCl,24,TColor::kRed);
 463   SetStyles(hAllProtonsOverLowPtProtonsDataCl,20,TColor::kRed);
 464
 465   //cerr<<"349"<<endl;
 466
 467   Rescale(hAntiNeutronsOverLowPtAntiProtons,2);
 468   Rescale(hNeutronsOverLowPtProtons,2);
 469   Rescale(hAntiProtonsOverLowPtAntiProtons,2);
 470   Rescale(hProtonsOverLowPtProtons,2);
 471   Rescale(hAllNeutronsOverLowPtProtonsData,2);
 472   Rescale(hAllNeutronsOverLowPtProtons,2);
 473   Rescale(hAllProtonsOverLowPtProtons,2);
 474   Rescale(hAllProtonsOverLowPtProtonsData,2);
 475
 476   Rescale(hAntiNeutronsOverLowPtAntiProtonsCl,2);
 477   Rescale(hNeutronsOverLowPtProtonsCl,2);
 478   Rescale(hAntiProtonsOverLowPtAntiProtonsCl,2);
 479   Rescale(hProtonsOverLowPtProtonsCl,2);
 480   Rescale(hAllNeutronsOverLowPtProtonsDataCl,2);
 481   Rescale(hAllNeutronsOverLowPtProtonsCl,2);
 482   Rescale(hAllProtonsOverLowPtProtonsCl,2);
 483   Rescale(hAllProtonsOverLowPtProtonsDataCl,2);
 484
 485
 486   //SetStyles(hAllProtonSim,28,1);
 487   //SetStyles(hAntiProtonData,22,TColor::kRed);
 488   //SetStyles(hProtonData,26,TColor::kRed);
 489   //SetStyles(hAllProtonData,33,1);
 490   hAntiNeutronsOverLowPtAntiProtons->SetMinimum(0.0);
 491   hAntiNeutronsOverLowPtAntiProtons->SetMaximum(8.0);
 492   hAntiNeutronsOverLowPtAntiProtons->GetXaxis()->SetTitle("N_{ch}");
 493   hAntiNeutronsOverLowPtAntiProtons->GetYaxis()->SetTitle("Ratio of energy deposited in sim");
 494   hAntiNeutronsOverLowPtAntiProtons->Draw("");
 495   hNeutronsOverLowPtProtons->Draw("same");
 496   hAntiProtonsOverLowPtAntiProtons->Draw("same");
 497   hProtonsOverLowPtProtons->Draw("same");
 498   hAllNeutronsOverLowPtProtons->Draw("same");
 499   hAllProtonsOverLowPtProtons->Draw("same");
 500   TLegend *legRatio2 = new TLegend(0.157718,0.709677,0.278523,0.935484);
 501   legRatio2->SetFillStyle(0);
 502   legRatio2->SetFillColor(0);
 503   legRatio2->SetBorderSize(0);
 504   legRatio2->SetTextSize(0.03);
 505   legRatio2->SetTextSize(0.038682);
 506   legRatio2->AddEntry(hAntiNeutronsOverLowPtAntiProtons,"#bar{n}/PID'd #bar{p}");
 507   legRatio2->AddEntry(hNeutronsOverLowPtProtons,"n/PID'd p");
 508   legRatio2->AddEntry(hAntiProtonsOverLowPtAntiProtons,"#bar{p}/PID'd #bar{p}");
 509   legRatio2->AddEntry(hProtonsOverLowPtProtons,"p/PID'd p");
 510   legRatio2->AddEntry(hAllNeutronsOverLowPtProtons,"(#bar{n}+n)/(PID'd #bar{p}+p)");
 511   legRatio2->AddEntry(hAllProtonsOverLowPtProtons,"(#bar{p}+p)/(PID'd #bar{p}+p)");
 512   legRatio2->Draw();
 513
 514   //cerr<<"387"<<endl;
 515   TString name4 = "/tmp/RatiosToLowPt"+detector+".eps";
 516   c4->SaveAs(name4.Data());
 517
 518   TCanvas *c5 = new TCanvas("c5","Neutron ratios used for error bounds",600,400);
 519   c5->SetTopMargin(0.02);
 520   c5->SetRightMargin(0.149329);
 521   c5->SetBorderSize(0);
 522   c5->SetFillColor(0);
 523   c5->SetFillColor(0);
 524   c5->SetBorderMode(0);
 525   c5->SetFrameFillColor(0);
 526   c5->SetFrameBorderMode(0);
 527
 528   TF1 *funcNominal = new TF1("funcNominal","([0]*exp(-x/[1])+[2]-[3]*x)*[4]",0,3700);
 529   funcNominal->SetParameter(0,3.38452e+00);
 530   funcNominal->SetParameter(1,1.32723e+02);
 531   funcNominal->SetParameter(2,2.51044e+00);
 532   funcNominal->SetParameter(3,4.78729e-04);
 533   funcNominal->SetParLimits(1,1,1e3);
 534   funcNominal->FixParameter(4,1);
 535   //cout<<"funcnominal"<<endl;
 536   funcNominal->SetLineColor(hAllNeutronsOverLowPtProtonsData->GetLineColor());
 537   hAllNeutronsOverLowPtProtonsData->Fit(funcNominal,"","",50,2700);
 538   TF1 *funcLow = funcNominal->Clone("funcLow");
 539   funcLow->FixParameter(4,1.15);
 540   funcLow->SetLineStyle(2);
 541   TF1 *funcHigh = funcNominal->Clone("funcHigh");
 542   funcHigh->FixParameter(4,0.85);
 543   funcHigh->SetLineStyle(2);
 544
 545
 546   hAllNeutronsOverLowPtProtons->SetMinimum(0.0);
 547   hAllNeutronsOverLowPtProtons->SetMaximum(8.0);
 548   hAllNeutronsOverLowPtProtons->GetXaxis()->SetTitle("N_{ch}");
 549   hAllNeutronsOverLowPtProtons->GetYaxis()->SetTitle("Ratio of energy deposited in sim");
 550   hAllNeutronsOverLowPtProtons->Draw();
 551   hAllNeutronsOverLowPtProtonsData->Draw("same");
 552   hAllProtonsOverLowPtProtons->Draw("same");
 553   hAllProtonsOverLowPtProtonsData->Draw("same");
 554   TLegend *legRatioForData = new TLegend(0.157718,0.709677,0.278523,0.935484);
 555   legRatioForData->SetFillStyle(0);
 556   legRatioForData->SetFillColor(0);
 557   legRatioForData->SetBorderSize(0);
 558   legRatioForData->SetTextSize(0.03);
 559   legRatioForData->SetTextSize(0.038682);
 560   legRatioForData->AddEntry(hAllNeutronsOverLowPtProtons,"(#bar{n}+n)/(PID'd #bar{p}+p in Sim)");
 561   legRatioForData->AddEntry(hAllNeutronsOverLowPtProtonsData,"(#bar{n}+n)/(PID'd #bar{p}+p in Data)");
 562   legRatioForData->AddEntry(hAllProtonsOverLowPtProtons,"(#bar{p}+p)/(PID'd #bar{p}+p in Sim)");
 563   legRatioForData->AddEntry(hAllProtonsOverLowPtProtonsData,"(#bar{p}+p)/(PID'd #bar{p}+p in Data)");
 564   legRatioForData->Draw();
 565   funcHigh->Draw("same");
 566   funcLow->Draw("same");
 567
 568   TString name5 = "/tmp/RatiosForErrors"+detector+".eps";
 569   c5->SaveAs(name5.Data());
 570
 571   TCanvas *c5a = new TCanvas("c5a","Neutron ratios used for error bounds",600,400);
 572   c5a->SetTopMargin(0.02);
 573   c5a->SetRightMargin(0.149329);
 574   c5a->SetBorderSize(0);
 575   c5a->SetFillColor(0);
 576   c5a->SetFillColor(0);
 577   c5a->SetBorderMode(0);
 578   c5a->SetFrameFillColor(0);
 579   c5a->SetFrameBorderMode(0);
 580
 581   TF1 *funcANominal = new TF1("funcANominal","([0]*exp(-x/[1])+[2]-[3]*x)*[4]",0,3700);
 582   funcANominal->SetParameter(0,3.09123e+00);
 583   funcANominal->SetParameter(1,2.69186e+01);
 584   funcANominal->SetParameter(2,3.07388e+00);
 585   funcANominal->SetParameter(3,2.88916e-03);
 586   funcANominal->SetParLimits(1,0,1e3);
 587   //funcANominal->SetParLimits(2,0,3);
 588   //funcANominal->SetParLimits(3,1e-3,1e3);
 589   funcANominal->SetParLimits(0,1,1e2);
 590   funcANominal->FixParameter(4,1);
 591   //cout<<"funcAnominal"<<endl;
 592   funcANominal->SetLineColor(hAllNeutronsOverLowPtProtonsDataCl->GetLineColor());
 593   Float_t fitlim = 280;
 594   if(isPhos) fitlim = 160;
 595   hAllNeutronsOverLowPtProtonsDataCl->Fit(funcANominal,"","",20,fitlim);
 596   TF1 *funcALow = funcANominal->Clone("funcALow");
 597   funcALow->FixParameter(4,1.15);
 598   funcALow->SetLineStyle(2);
 599   TF1 *funcAHigh = funcANominal->Clone("funcAHigh");
 600   funcAHigh->FixParameter(4,0.85);
 601   funcAHigh->SetLineStyle(2);
 602
 603
 604   hAllNeutronsOverLowPtProtonsCl->SetMinimum(0.0);
 605   hAllNeutronsOverLowPtProtonsCl->SetMaximum(8.0);
 606   hAllNeutronsOverLowPtProtonsCl->GetXaxis()->SetTitle("N_{cl}");
 607   hAllNeutronsOverLowPtProtonsCl->GetYaxis()->SetTitle("Ratio of energy deposited in sim");
 608   hAllNeutronsOverLowPtProtonsCl->GetXaxis()->SetRange(1,hAllNeutronsOverLowPtProtonsCl->FindBin(fitlim));
 609   hAllNeutronsOverLowPtProtonsCl->Draw();
 610   hAllNeutronsOverLowPtProtonsDataCl->Draw("same");
 611   hAllProtonsOverLowPtProtonsCl->Draw("same");
 612   hAllProtonsOverLowPtProtonsDataCl->Draw("same");
 613   TLegend *legRatioForDataCl = new TLegend(0.157718,0.709677,0.278523,0.935484);
 614   legRatioForDataCl->SetFillStyle(0);
 615   legRatioForDataCl->SetFillColor(0);
 616   legRatioForDataCl->SetBorderSize(0);
 617   legRatioForDataCl->SetTextSize(0.03);
 618   legRatioForDataCl->SetTextSize(0.038682);
 619   legRatioForDataCl->AddEntry(hAllNeutronsOverLowPtProtonsCl,"(#bar{n}+n)/(PID'd #bar{p}+p in Sim)");
 620   legRatioForDataCl->AddEntry(hAllNeutronsOverLowPtProtonsDataCl,"(#bar{n}+n)/(PID'd #bar{p}+p in Data)");
 621   legRatioForDataCl->AddEntry(hAllProtonsOverLowPtProtonsCl,"(#bar{p}+p)/(PID'd #bar{p}+p in Sim)");
 622   legRatioForDataCl->AddEntry(hAllProtonsOverLowPtProtonsDataCl,"(#bar{p}+p)/(PID'd #bar{p}+p in Data)");
 623   legRatioForDataCl->Draw();
 624   funcAHigh->Draw("same");
 625   funcALow->Draw("same");
 626
 627   TString name5a = "/tmp/RatiosForErrorsCl"+detector+".eps";
 628   c5a->SaveAs(name5a.Data());
 629
 630
 631
 632 //     TCanvas *c5b = new TCanvas("c5b","Neutron ratios used for error bounds",600,400);
 633 //     c5b->SetTopMargin(0.02);
 634 //     c5b->SetRightMargin(0.149329);
 635 //     c5b->SetBorderSize(0);
 636 //     c5b->SetFillColor(0);
 637 //     c5b->SetFillColor(0);
 638 //     c5b->SetBorderMode(0);
 639 //     c5b->SetFrameFillColor(0);
 640 //     c5b->SetFrameBorderMode(0);
 641 //     fHistNeutronsDepositedVsNclProf->Draw();
 642 //     fHistAntiNeutronsDepositedVsNclProf->Draw("same");
 643 //     fHistPIDProtonsTrackMatchedDepositedVsNclProf->Draw("same");
 644 //     fHistPIDAntiProtonsTrackMatchedDepositedVsNclProf->Draw("same");
 645 //     //fHistPIDProtonsTrackMatchedDepositedVsNclProf->Draw("same");
 646 //     //fHistPIDAntiProtonsTrackMatchedDepositedVsNclProf->Draw("same");
 647
 648 //     SetStyles(fHistNeutronsDepositedVsNclProf,20, TColor::kBlue);
 649 //     SetStyles(fHistAntiNeutronsDepositedVsNclProf,24, TColor::kBlue);
 650 //     SetStyles(fHistPIDProtonsTrackMatchedDepositedVsNclProf,21,TColor::kRed);
 651 //     SetStyles(fHistPIDAntiProtonsTrackMatchedDepositedVsNclProf,25,TColor::kRed);
 652 //     PrintInfo(fHistNeutronsDepositedVsNclProf);
 653 //     PrintInfo(fHistAntiNeutronsDepositedVsNclProf);
 654 //     PrintInfo(fHistPIDProtonsTrackMatchedDepositedVsNclProf);
 655 //     PrintInfo(fHistPIDAntiProtonsTrackMatchedDepositedVsNclProf);
 656
 657     //cerr<<"442"<<endl;
 658
 659     TObjArray trackmultiplicity(20);
 660     TObjArray trackmultiplicityData(20);
 661     TObjArray clustermultiplicity(20);
 662     TObjArray clustermultiplicityData(20);
 663     int nbinsChMult = fHistCentVsNchVsNcl->GetYaxis()->GetNbins();
 664     int nbinsClMult = fHistCentVsNchVsNcl->GetZaxis()->GetNbins();
 665     fHistCentVsNchVsNcl->GetXaxis()->SetTitle("Cent Bin");
 666     fHistCentVsNchVsNcl->GetYaxis()->SetTitle("N_{tr}");
 667     fHistCentVsNchVsNcl->GetZaxis()->SetTitle("N_{cl}");
 668     for(int cb=0;cb<20;cb++){
 669       //x axis is centrality
 670       //y axis is charged track multiplicity
 671       //z axis is cluster multiplicity
 672       fHistCentVsNchVsNcl->GetXaxis()->SetRange(cb+1,cb+1);
 673       trackmultiplicity[cb] = fHistCentVsNchVsNcl->Project3D("y");
 674       SetStyles((TH1*)trackmultiplicity[cb],markers[cb],colors[cb],Form("tr%i",cb));
 675       clustermultiplicity[cb] = fHistCentVsNchVsNcl->Project3D("z");
 676       SetStyles((TH1*)clustermultiplicity[cb],markers[cb],colors[cb],Form("cl%i",cb));
 677       fHistCentVsNchVsNclData->GetXaxis()->SetRange(cb+1,cb+1);
 678       trackmultiplicityData[cb] = fHistCentVsNchVsNclData->Project3D("y");
 679       SetStyles((TH1*)trackmultiplicityData[cb],markers[cb],colors[cb],Form("tr%i",cb));
 680       clustermultiplicityData[cb] = fHistCentVsNchVsNclData->Project3D("z");
 681       SetStyles((TH1*)clustermultiplicityData[cb],markers[cb],colors[cb],Form("cl%i",cb));
 682     }
 683
 684     int currentShortCB = 0;
 685     //cerr<<"464"<<endl;
 686     Float_t neventsShortData[] = {0,0,0,0,0,0,0,0,0,0,0};
 687     Float_t neventsShortSim[] = {0,0,0,0,0,0,0,0,0,0,0};
 688     Float_t neventsShortDataCl[] = {0,0,0,0,0,0,0,0,0,0,0};
 689     Float_t neventsShortSimCl[] = {0,0,0,0,0,0,0,0,0,0,0};
 690     for(int cb=0;cb<19;cb++){
 691       //cout<<"cb "<<cb<<" current short cb "<<currentShortCB<<endl;
 692       cout<<"cb "<<cb<<" mean mult "<< ((TH1*)trackmultiplicity[cb])->GetMean()<<endl;
 693       int nbins = ((TH1*)trackmultiplicity[cb])->GetNbinsX();
 694       Float_t nevents = 0.0;
 695       for(int binn = 1;binn<=nbins;binn++){
 696         float neventsBinn = ((TH1*)trackmultiplicityData[cb])->GetBinContent(binn);
 697         float mult= ((TH1*)trackmultiplicityData[cb])->GetBinCenter(binn);
 698         //mean et deposited in this event class =  (et of n+nbar / measured et of p+pbar) * (measured et of p+pbar)
 699         float meanetppbar = fHistPIDAntiProtonsTrackMatchedDepositedVsNchDataProf->GetBinContent(fHistPIDAntiProtonsTrackMatchedDepositedVsNchDataProf->FindBin(mult)) +  fHistPIDProtonsTrackMatchedDepositedVsNchDataProf->GetBinContent(fHistPIDProtonsTrackMatchedDepositedVsNchDataProf->FindBin(mult));
 700         float meanet = funcNominal->Eval(mult) * meanetppbar;
 701         nNeutronsData[cb] += neventsBinn*meanet;
 702         nevents +=neventsBinn;
 703         nNeutronsShortData[currentShortCB] += neventsBinn*meanet;
 704         neventsShortData[currentShortCB] +=neventsBinn;
 705       }
 706       if(nevents>0) nNeutronsData[cb] = nNeutronsData[cb]/nevents;
 707       nevents = 0.0;
 708       for(int binn = 1;binn<=nbins;binn++){
 709         float neventsBinn = ((TH1*)trackmultiplicity[cb])->GetBinContent(binn);
 710         float mult= ((TH1*)trackmultiplicity[cb])->GetBinCenter(binn);
 711         //mean et deposited in this event class =  (et of n+nbar / measured et of p+pbar) * (measured et of p+pbar)
 712         float meanetppbar = fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf->GetBinContent(fHistPIDAntiProtonsTrackMatchedDepositedVsNchProf->FindBin(mult)) +  fHistPIDProtonsTrackMatchedDepositedVsNchProf->GetBinContent(fHistPIDProtonsTrackMatchedDepositedVsNchProf->FindBin(mult));
 713         float meanet = funcNominal->Eval(mult) * meanetppbar;
 714         nNeutronsSim[cb] += neventsBinn*meanet;
 715         nevents +=neventsBinn;
 716         nNeutronsShortSim[currentShortCB] += neventsBinn*meanet;
 717         neventsShortSim[currentShortCB] +=neventsBinn;
 718       }
 719       if(nevents>0) nNeutronsSim[cb] = nNeutronsSim[cb]/nevents;
 720
 721
 722       nbins = ((TH1*)clustermultiplicity[cb])->GetNbinsX();
 723       nevents = 0.0;
 724       for(int binn = 1;binn<=nbins;binn++){
 725         float neventsBinn = ((TH1*)clustermultiplicityData[cb])->GetBinContent(binn);
 726         float mult= ((TH1*)clustermultiplicityData[cb])->GetBinCenter(binn);
 727         //mean et deposited in this event class =  (et of n+nbar / measured et of p+pbar) * (measured et of p+pbar)
 728         float meanetppbar = fHistPIDAntiProtonsTrackMatchedDepositedVsNclDataProf->GetBinContent(fHistPIDAntiProtonsTrackMatchedDepositedVsNclDataProf->FindBin(mult)) +  fHistPIDProtonsTrackMatchedDepositedVsNclDataProf->GetBinContent(fHistPIDProtonsTrackMatchedDepositedVsNclDataProf->FindBin(mult));
 729         float meanet = funcANominal->Eval(mult) * meanetppbar;
 730         nNeutronsDataCl[cb] += neventsBinn*meanet;
 731         nevents +=neventsBinn;
 732         //cout<<"et "<< neventsBinn<<"*"<<meanet<<" events "<<neventsBinn<<endl;
 733         nNeutronsShortDataCl[currentShortCB] += neventsBinn*meanet;
 734         neventsShortDataCl[currentShortCB] +=neventsBinn;
 735       }
 736       if(nevents>0) nNeutronsDataCl[cb] = nNeutronsDataCl[cb]/nevents;
 737       nevents = 0.0;
 738       for(int binn = 1;binn<=nbins;binn++){
 739         float neventsBinn = ((TH1*)clustermultiplicity[cb])->GetBinContent(binn);
 740         float mult= ((TH1*)clustermultiplicity[cb])->GetBinCenter(binn);
 741         //mean et deposited in this event class =  (et of n+nbar / measured et of p+pbar) * (measured et of p+pbar)
 742         float meanetppbar = fHistPIDAntiProtonsTrackMatchedDepositedVsNclProf->GetBinContent(fHistPIDAntiProtonsTrackMatchedDepositedVsNclProf->FindBin(mult)) +  fHistPIDProtonsTrackMatchedDepositedVsNclProf->GetBinContent(fHistPIDProtonsTrackMatchedDepositedVsNclProf->FindBin(mult));
 743         float meanet = funcANominal->Eval(mult) * meanetppbar;
 744         nNeutronsSimCl[cb] += neventsBinn*meanet;
 745         nevents +=neventsBinn;
 746         nNeutronsShortSimCl[currentShortCB] += neventsBinn*meanet;
 747         neventsShortSimCl[currentShortCB] +=neventsBinn;
 748       }
 749       if(nevents>0) nNeutronsSimCl[cb] = nNeutronsSimCl[cb]/nevents;
 750
 751
 752       //cout<<"cb "<<cb;
 753       //cout<<" data "<<  nNeutronsData[cb];// <<" +/- "<< 0.15*nNeutronsData[cb];
 754       //cout<<" sim "<<  nNeutronsSim[cb];// <<" +/- "<< 0.15*nNeutronsSim[cb]<<endl;
 755       //cout<<" data cl "<<  nNeutronsDataCl[cb];// <<" +/- "<< 0.15*nNeutronsDataCl[cb];
 756       //cout<<" sim cl "<<  nNeutronsSimCl[cb];// <<" +/- "<< 0.15*nNeutronsSimCl[cb]<<endl;
 757       //cout<<" w/err data ";
 758       float val = (nNeutronsData[cb]+nNeutronsDataCl[cb])/2.0;
 759       float err = TMath::Abs(nNeutronsData[cb]-nNeutronsDataCl[cb])/2.0;
 760       //cout<< val<<" +/- "<<err;
 761       //if(val>0)cout<<" ("<<err/val<<")";
 762       float valsim = (nNeutronsSim[cb]+nNeutronsSimCl[cb])/2.0;
 763       float errsim = TMath::Abs(nNeutronsSim[cb]-nNeutronsSimCl[cb])/2.0;
 764       //cout<<" sim "<< valsim <<" +/- "<<errsim;
 765       //if(valsim>0)cout<<" ("<<errsim/valsim<<")";
 766       //cout<<endl;
 767       nNeutronsData[cb] = val;
 768       nNeutronsDataErr[cb] = err;
 769       myfile<<Form("%2.3f %2.3f",nNeutronsData[cb],nNeutronsDataErr[cb])<<endl;
 770       if(cb<2 || cb%2==1){//normalize
 771         if(neventsShortSimCl[currentShortCB]>0) nNeutronsShortSimCl[currentShortCB] = nNeutronsShortSimCl[currentShortCB]/neventsShortSimCl[currentShortCB];
 772         if(neventsShortDataCl[currentShortCB]>0) nNeutronsShortDataCl[currentShortCB] = nNeutronsShortDataCl[currentShortCB]/neventsShortDataCl[currentShortCB];
 773         if(neventsShortSim[currentShortCB]>0) nNeutronsShortSim[currentShortCB] = nNeutronsShortSim[currentShortCB]/neventsShortSim[currentShortCB];
 774         if(neventsShortData[currentShortCB]>0) nNeutronsShortData[currentShortCB] = nNeutronsShortData[currentShortCB]/neventsShortData[currentShortCB];
 775         //cout<<"cbShort "<<currentShortCB<<" data "<<  nNeutronsShortData[currentShortCB] <<" +/- "<< 0.15*nNeutronsShortData[currentShortCB];
 776         //cout<<" sim "<<  nNeutronsShortSim[currentShortCB] <<" +/- "<< 0.15*nNeutronsShortSim[currentShortCB];//<<endl;
 777
 778         val = (nNeutronsShortData[currentShortCB]+nNeutronsShortDataCl[currentShortCB])/2.0;
 779         err = TMath::Abs(nNeutronsShortData[currentShortCB]-nNeutronsShortDataCl[currentShortCB])/2.0;
 780         //cout<<" val "<<val<<" err "<<err<<endl;
 781         //cout<<nNeutronsShortData[currentShortCB]<<" "<<nNeutronsShortDataCl[currentShortCB]<<endl<<endl;
 782         nNeutronsShortData[currentShortCB] = val;
 783         nNeutronsShortDataErr[currentShortCB] = err;
 784
 785
 786         myfile2<<Form("%2.3f %2.3f",nNeutronsShortData[currentShortCB],0.15*nNeutronsShortData[currentShortCB])<<endl;
 787         if(cb<2) currentShortCB++;
 788       }
 789       if(cb>2 && cb%2==1){//increment the counter
 790         //cout<<"Here!"<<endl;
 791         currentShortCB++;
 792       }
 793
 794     }
 795     myfile.close();
 796     myfile2.close();
 797     WriteLatex();
 798 }
 799
 800
 801 void WriteLatex(){
 802   TString detector = "Emcal";
 803     string inline;
 804     //NeutronsEmcal.dat
 805     TString neutronInfileName = "Neutrons"+detector+".dat";
 806     ifstream myneutronfile3 (neutronInfileName.Data());
 807     Float_t value = 0;
 808     Float_t error = 0;
 809     Int_t i=0;
 810     if (myneutronfile3.is_open()){
 811       while ( myneutronfile3.good() )
 812         {
 813           getline (myneutronfile3,inline);
 814           istringstream tmp(inline);
 815           tmp >> value;
 816           tmp >> error;
 817           if(i<20){
 818             neutronCorrEmcal[i] = value;
 819             neutronErrorEmcal[i] = error;
 820           }
 821           i++;
 822         }
 823         myneutronfile3.close();
 824     }
 825
 826     detector = "Phos";
 827     neutronInfileName = "Neutrons"+detector+".dat";
 828     ifstream myneutronfile4 (neutronInfileName.Data());
 829     Float_t value = 0;
 830     Float_t error = 0;
 831     Int_t i=0;
 832     if (myneutronfile4.is_open()){
 833       while ( myneutronfile4.good() )
 834         {
 835           getline (myneutronfile4,inline);
 836           istringstream tmp(inline);
 837           tmp >> value;
 838           tmp >> error;
 839           if(i<20){
 840             neutronCorrPhos[i] = value;
 841             neutronErrorPhos[i] = error;
 842           }
 843           i++;
 844         }
 845         myneutronfile4.close();
 846     }
 847     ofstream myfile3;
 848     myfile3.open ("Protons.tex");
 849     for(int i=0;i<20;i++){
 850       TString line = Form("%i-%i & %2.3f $\\pm$ %2.3f & %2.3f $\\pm$ %2.3f \\\\",i*5,(i+1)*5,neutronCorrPhos[i],neutronErrorPhos[i],neutronCorrEmcal[i],neutronErrorEmcal[i]);
 851       myfile3<<line.Data()<<endl;
 852
 853     }
 854
 855     myfile3.close();
 856
 857 }
 858
 859
 860 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name)
 861 {
 862     if(!numerator){
 863       cerr<<"Error:  numerator does not exist!"<<endl;
 864       return NULL;
 865     }
 866     if(!denominator){
 867       cerr<<"Error:  denominator does not exist!"<<endl;
 868       return NULL;
 869     }
 870     TH1* result = (TH1*)numerator->Clone(name);
 871     Int_t nbins = numerator->GetNbinsX();
 872     for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
 873       Double_t numeratorVal = numerator->GetBinContent(ibin);
 874       Double_t denominatorVal = denominator->GetBinContent(ibin);
 875       // Check if the errors are right or the thing is scaled
 876       Double_t numeratorValErr = numerator->GetBinError(ibin);
 877       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
 878         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
 879         numeratorVal /= rescale;
 880       }
 881       Double_t denominatorValErr = denominator->GetBinError(ibin);
 882       if (!(denominatorValErr==0. || denominatorVal==0. )) {
 883         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
 884         denominatorVal /= rescale;
 885       }
 886       Double_t quotient = 0.;
 887       if (denominatorVal!=0.) {
 888         quotient = numeratorVal/denominatorVal;
 889       }
 890       Double_t quotientErr=0;
 891       quotientErr = TMath::Sqrt(
 892                                 (numeratorVal+1.0)/(denominatorVal+2.0)*
 893                                 ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
 894       result->SetBinContent(ibin,quotient);
 895       result->SetBinError(ibin,quotientErr);
 896       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
 897     }
 898     return result;
 899 }