PWGLF/totEt/macros/emEt/kaonCorr/SecondaryK0SJacSys.C

   1 #include <iostream>
   2 #include <fstream>
   3 #include "TMath.h"
   4
   5 using namespace std;
   6
   7 char myjunk[200];
   8
   9 void SetStyles(TH1 *histo, char *ytitle, char *xtitle, int color, int marker);
  10 void SetCanvasStyle(TCanvas *c);
  11 TLegend *GetLegend(float x1, float y1, float x2, float y2);
  12 Bool_t doLevy = kFALSE;//for debugging
  13 Bool_t doBlast = kTRUE;//for debugging
  14 void Draw(TH1 *data, TF1 *fLevy, TF1 *fBlast,int centbin,char *pid,char *name);
  15 void PrintLatex(int etCutNum, char *det,Bool_t effCorr);
  16 void PrintArrays();
  17
  18 // Centrality dependent factors
  19 int nbins = 10;  // number of centrality bins
  20 Double_t etk0[2][10] = {{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}};
  21 Double_t etBlastk0[2][10] = {{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}};
  22 Double_t etBlastk02D[2][10] = {{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}};
  23 //et cuts for kaon energies
  24 //0.00, 0.050, 0.100, 0.150, 0.200
  25 //0.250,0.300, 0.350, 0.400, 0.450
  26 //0.500
  27 //counting is integer and starts are 1
  28 //                        1     2      3      4      5       6PHOS 7EMCal 8      9      10      11Alt
  29 Double_t etCutOffs[11] = {0.00, 0.050, 0.100, 0.150, 0.200,  0.250,0.300, 0.350, 0.400, 0.450,  0.500};
  30
  31 Double_t kaonDeposits[2][10] = {{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}};
  32
  33 //note:  Not at all sensitive to which Jacobian is used, which also means the pT used doesn't really matter
  34 void SecondaryK0SJacSys(char *inputfilename = "../workingdir/rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root", char *det = "EMCal", int whichjacobian = 1, int etCutNum = 7, Bool_t longRun = kFALSE, Bool_t effCorr = kFALSE){//Kaon collection: 0=K0S, -1=K-, +1=K+
  35   Int_t kaonSelection = 1;
  36   gStyle->SetOptTitle(0);
  37   gStyle->SetOptStat(0);
  38   gStyle->SetOptFit(0);
  39   gStyle->SetLineWidth(3);
  40   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  41   gROOT->ProcessLine(".L AliBWFunc.cxx+g");
  42   gROOT->ProcessLine(".L AliBWTools.cxx+g");
  43   gSystem->Load("libTree.so");
  44   gSystem->Load("libVMC.so");
  45   gSystem->Load("libMinuit.so");
  46   gSystem->Load("libSTEERBase.so");
  47   gSystem->Load("libESD.so");
  48   gSystem->Load("libAOD.so");
  49   gSystem->Load("libANALYSIS.so");
  50   gSystem->Load("libANALYSISalice.so");
  51   gSystem->Load("libCORRFW.so");
  52   gSystem->Load("libPWGLFspectra.so");
  53   gSystem->Load("libPWGTools.so");
  54   gROOT->LoadMacro("macros/FitParticle.C+");
  55   gROOT->LoadMacro("macros/GetLevidEtdptTimesPt.C");
  56   //inputting kaon errors from http://arxiv.org/abs/1303.0737
  57   //using average of K+ and K-
  58   Float_t kaonYield[10] = {109,90.5,68,46,30,18.2,10.2,5.1,2.3,0.855};
  59   Float_t kaonStatErr[10] = {0.3,0.2,0.1,0.1,0.1, 0.06,0.04,0.03,0.02,0.01};
  60   Float_t kaonSysErr[10] = {9,7,5,4,2, 1.5,0.8,0.4,0.2,0.09};
  61   Float_t kaonTotErr[10] = {0,0,0,0,0, 0,0,0,0,0};
  62   Float_t kaonFracErr[10] = {0,0,0,0,0, 0,0,0,0,0};
  63   for(int i=0;i<10;i++){
  64     kaonTotErr[i] = TMath::Sqrt(kaonSysErr[i]*kaonSysErr[i]+kaonStatErr[i]*kaonStatErr[i]);
  65     kaonFracErr[i] = kaonTotErr[i]/kaonYield[i];
  66     //cout<<"bin "<<i<<" err "<<kaonFracErr[i]<<endl;
  67   }
  68
  69   //centrality bin #  (want to loop over 5 bins)
  70   int centbin = 0;  // can change later
  71   TString detector = det;
  72   float etarange = 0.7;
  73   if(detector.Contains("EMC")){
  74     etarange = 0.7;
  75   }
  76   else{etarange = 0.1;}
  77
  78   //==================READ IN DATA==========================
  79
  80   TFile *chargedkaonfile = new TFile("rootFiles/SPECTRA_COMB_20120709.root");
  81   gROOT->LoadMacro("macros/k0sFinal.C");                // load data spectra
  82   TClonesArray histosk0 = GetK0S();                     // get K0S histogram
  83
  84   TCanvas *cChK = new TCanvas("cChK","ChargedKaons",700,500);
  85   //cChK->Divide(4,3);
  86   TH1D *histoskch[] = {NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL, NULL,NULL};
  87   if(kaonSelection==1){
  88     for(int i=0;i<=9;i++){
  89       TH1D *histoTemp = (TH1D*) chargedkaonfile->Get(Form("cent%i_kaon_plus",i));
  90       histoskch[i] = (TH1D*)histoTemp->Clone("cent%i_kaon_plus_clone");
  91       if(i==0){
  92         ((TH1D*)histoskch[i])->Draw();
  93       }
  94       else{
  95         ((TH1D*)histoskch[i])->Draw("same");
  96       }
  97       //}
  98     }
  99   }
 100   else{
 101     for(int i=0;i<=9;i++){
 102       TH1D *histoTemp = (TH1D*) chargedkaonfile->Get(Form("cent%i_kaon_minus",i));
 103       histoskch[i] = (TH1D*)histoTemp->Clone("cent%i_kaon_minus_clone");
 104       if(i==0){
 105         ((TH1D*)histoskch[i])->Draw();
 106       }
 107       else{
 108         ((TH1D*)histoskch[i])->Draw("same");
 109       }
 110       //}
 111     }
 112   }
 113   TCanvas *cEmpty = new TCanvas("cEmpty","Empty",700,500);
 114
 115   TString particleName = "K0";
 116   TString particleLatex = "K^{0}_S";
 117   if(kaonSelection==1){
 118     particleName = "K+";
 119     particleLatex = "K^{+}";
 120   }
 121   if(kaonSelection==-1){
 122     particleName = "K-";
 123     particleLatex = "K^{-}";
 124   }
 125
 126   //==================FIT DATA==========================
 127   TF1 *funck0[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};
 128   TF1 *funcBlastk0[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};     // Blastwave fit to K0 pT?
 129   TF1 *funcetk0[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};
 130   TF1 *funcetBlastk0[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};   // Blastwave fit to K0 et
 131   TF1 *funcetBlastk02D[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};
 132   TF1 *funcBlastk0Jac[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};
 133
 134   TH1 *histoSysUp[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};  //K0S Clone histogram
 135   TF1 *funcBlastk0Up[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};
 136   TH1 *histoSysLow[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};  //K0S Clone histogram
 137   TF1 *funcBlastk0Low[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};
 138
 139   TH1 *histoSys[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};  //K0S Clone histogram
 140
 141
 142   Float_t TInit = 0.1326;    // initialize T parameter of Blastwave fit
 143   Float_t nInit = 6.6857;    // initialize n parameter of Blastwave fit
 144   Float_t normInit = 20;     // set normalization parameter of Blastwave fit
 145   AliBWTools *tools = new AliBWTools();
 146   Double_t yield;            // yield
 147   Double_t yieldError;       // yield error
 148   Double_t partialYields[] = {0,0,0};
 149   Double_t partialYieldErrors[] = {0,0,0};
 150   Double_t Ameson = 0;
 151   Double_t Abaryon = -1;
 152   Double_t Aantibaryon = 1;
 153
 154   //Systematic erros
 155
 156   TF1 *funcBlastk0Par1[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};
 157   TF1 *funcBlastk0Par2[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};
 158   Double_t Par[10] = {0,0,0,0,0, 0,0,0,0,0};
 159   Double_t ErrPar[10] = {0,0,0,0,0, 0,0,0,0,0};
 160
 161   // Genarate y
 162
 163   TRandom* randy = new TRandom();
 164
 165   // integral of fit to data
 166   float fitINTEGRALdata[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};
 167
 168   if(kaonSelection!=0){//charged kaons
 169     nbins = 10;
 170   }
 171
 172   // loop over centrality bins 0->5 and genarate histos shift up and down to Systematic errors
 173   for(int i=0;i<nbins;i++){
 174     //for(int i=centbin;i<centbin+1;i++){
 175     if(kaonSelection==0){
 176       histoSys[i] = (TH1F*)histosk0[i]->Clone(Form("Sys%i",i));
 177       histoSysLow[i] = (TH1F*)histosk0[i]->Clone(Form("SysLow%i",i));
 178       histoSysUp[i] = (TH1F*)histosk0[i]->Clone(Form("SysUp%i",i));
 179     }
 180     else{
 181       histoSys[i] = (TH1D*)histoskch[i]->Clone(Form("Sys%i",i));
 182       histoSysLow[i] = (TH1D*)histoskch[i]->Clone(Form("SysLow%i",i));
 183       histoSysUp[i] = (TH1D*)histoskch[i]->Clone(Form("SysUp%i",i));
 184     }
 185     float NbinsSys =  histoSysUp[i]->GetXaxis()->GetNbins();
 186
 187     for(int x = 1; x<=NbinsSys; x++){
 188
 189       float BinContent = histoSysUp[i]->GetBinContent(x);
 190       float BinError =  histoSysUp[i]->GetBinError(x);
 191       histoSysUp[i]->SetBinContent(x, BinContent+BinError);
 192       histoSysLow[i]->SetBinContent(x, BinContent-BinError);
 193     }
 194
 195     funcBlastk0Up[i] = FitParticle((TH1*)histoSysUp[i],particleName.Data(),-1,-1,-1,0,0,3,-1,kFitBlastWave);
 196     funcBlastk0Low[i] = FitParticle((TH1*)histoSysLow[i],particleName.Data(),-1,-1,-1,0,0,3,-1,kFitBlastWave);
 197     funcBlastk0[i] = FitParticle((TH1*)histoSys[i],particleName.Data(),-1,-1,-1,0,0,3,-1,kFitBlastWave);
 198
 199     if(kaonSelection==0){
 200       funck0[i] = FitParticle((TH1F*)histosk0[i],particleName.Data(),TInit,normInit,nInit);
 201       if(!histosk0[i]) cerr<<"histosk0[i] does not exist!"<<endl;
 202       if(!funck0[i]) cerr<<"funck0[i] does not exist!"<<endl;
 203       Draw((TH1F*)histosk0[i],funck0[i],funcBlastk0[i],i,particleLatex.Data(),particleName.Data());
 204     }
 205     else{
 206       if(!histoskch[i]) cerr<<"Warning!  Histogram does not exist!"<<endl;
 207       funck0[i] = FitParticle((TH1D*)histoskch[i],particleName.Data());
 208       cerr<<"174"<<endl;
 209       Draw(histoskch[i],funck0[i],funcBlastk0[i],i,particleLatex.Data(),particleName.Data());
 210     }
 211   }
 212
 213
 214   //==================THROWING RANDOM SPECTRA==========================
 215   // number of particles to put in histogram (will change later)
 216   int nParticles = 1e3; // for quick runs
 217   if(longRun) nParticles = 1e8;   //for statistics
 218
 219   // set parameters of setstyles
 220   int marker1 = 21; int marker2 = 24;
 221   int color1 = 2; int color2 = 4;
 222   int mixmarker = 22; int mixcolor = 3;
 223
 224   TFile *f = TFile::Open(inputfilename, "READ");
 225   TList *l = (TList*)f->Get("out1");
 226   TH3F *fHistK0EDepositsVsPtInAcceptance;
 227   TH3F *fHistK0EDepositsVsPtOutOfAcceptance;
 228   if(effCorr){
 229     fHistK0EDepositsVsPtInAcceptance =(TH3F*) l->FindObject("fHistK0EDepositsVsPtInAcceptance");
 230     fHistK0EDepositsVsPtOutOfAcceptance =(TH3F*) l->FindObject("fHistK0EDepositsVsPtOutOfAcceptance");
 231   }
 232   else{
 233     fHistK0EDepositsVsPtInAcceptance =(TH3F*) l->FindObject("fHistK0EGammaVsPtInAcceptance");
 234     fHistK0EDepositsVsPtOutOfAcceptance =(TH3F*) l->FindObject("fHistK0EGammaVsPtOutOfAcceptance");
 235   }
 236   if(!fHistK0EDepositsVsPtOutOfAcceptance) cerr<<"Warning!  did not find fHistK0EDepositsVsPtOutOfAcceptance"<<endl;
 237   if(!fHistK0EDepositsVsPtInAcceptance) cerr<<"Warning!  did not find fHistK0EDepositsVsPtInAcceptance"<<endl;
 238   fHistK0EDepositsVsPtInAcceptance->GetZaxis()->SetRange(etCutNum,etCutNum);
 239   TH2D *my2DhistoK0S = (TH2D*) fHistK0EDepositsVsPtInAcceptance->Project3D("yx");
 240   fHistK0EDepositsVsPtOutOfAcceptance->GetZaxis()->SetRange(etCutNum,etCutNum);
 241   TH2D *my2DhistoK0SOutOfAcceptance = (TH2D*) fHistK0EDepositsVsPtOutOfAcceptance->Project3D("yx");
 242   TH1F *fHistSimKaonsInAcceptance = l->FindObject("fHistSimKaonsInAcceptance");
 243   TH1F *fHistSimKaonsOutOfAcceptance = l->FindObject("fHistSimKaonsOutOfAcceptance");
 244   TH1F *hRatioOutOfAccOverInAcc = fHistSimKaonsOutOfAcceptance->Clone("hRatioOutOfAccOverInAcc");
 245   hRatioOutOfAccOverInAcc->Divide(fHistSimKaonsInAcceptance);
 246
 247   //find number of bins of Caio's histogram
 248   int bincount = 0;
 249   bincount = my2DhistoK0S->GetXaxis()->GetNbins();
 250   cout<<"Number of bins in histo: "<<bincount<<endl;
 251
 252   // create array to hold histograms
 253   TClonesArray histoarrayK0S("TH1D",bincount);  //(nbins);
 254   TClonesArray histoarrayK0SOutOfAcceptance("TH1D",bincount);  //(nbins);
 255
 256   //==================PROJECTIONS FOR INPUT ET DEPOSITS==========================
 257   //loop over bins in x to define projections  (nbins)
 258   for(int i=0; i<bincount; i++){
 259
 260     // make projections of 2D histogram my2DhistoK0S
 261     histoarrayK0S[i] = (TH1D*)my2DhistoK0S->ProjectionY(Form("tmpK0S%i", i+1), i+1, i+1);
 262     histoarrayK0SOutOfAcceptance[i] = (TH1D*)my2DhistoK0SOutOfAcceptance->ProjectionY(Form("tmpK0SOutOfAcc%i", i+1), i+1, i+1);
 263   }
 264
 265
 266   for(int j=0;j<nbins;j++){
 267     centbin = j;
 268     cout<<"Working on centrality bin "<<centbin<<endl;
 269
 270
 271     // creates histogram of spectra from thrown K0S
 272     TH1F *histoSpectrumK0S = new TH1F("histoSpectrumK0S","p_{T} spectrum of K^{0}_{S}",100,0.0,5.0);
 273     SetStyles(histoSpectrumK0S,"dN/dp_{T}","p_{T} (GeV)", color1, marker1);
 274     TH1F *histoETspectrumK0S = new TH1F("histoETspectrumK0S","transverse energy scaled of K^{0}_{S}",100,0.0,5.0);
 275     SetStyles(histoETspectrumK0S,"dN/dE_{T} calculated","E_{T} (GeV)", color2, marker2);
 276     TH1F *histoETCspectrumK0S = new TH1F("histoETCspectrumK0S","transverse energy scaled from data of K^{0}_{S}",100,0.0,5.0);
 277     SetStyles(histoETCspectrumK0S,"dN/dE_{T} from data","E_{T} (GeV)", mixcolor, mixmarker);
 278
 279     // initialize totals
 280     float totET = 0;               // total ET of K0S's calculated
 281     float totetcK0S = 0;           // total ET of K0S's from data
 282     float testTOT = 0;
 283     float testTOTpt = 0;
 284     float pK0S = 0;
 285
 286     // declare edges of first and last bin
 287     int mybin100 = my2DhistoK0S->GetXaxis()->FindBin(0.1001);
 288     int mybin0 = my2DhistoK0S->GetXaxis()->FindBin(0.0001);
 289
 290     //Jacobian correction
 291
 292     TH1F *histoSpectrumK0SJac;
 293     histoSpectrumK0SJac = (TH1F*)histoSys[centbin]->Clone("Jac");
 294     SetStyles(histoSpectrumK0SJac,"dN/dp_{T}","p_{T} (GeV)", color1, marker1);
 295
 296     TH1F *histoSpectrumK0SDivide = (TH1F*)histoSys[centbin]->Clone("JacDiv");
 297
 298
 299 //     TH1F *histoeta = new TH1F("histoeta","",100,-1,1);
 300 //     SetStyles(histoeta,"N","#eta", color1, marker1);
 301
 302     //Jacobian Correction
 303
 304     float NbinsJac = histoSys[centbin]->GetXaxis()->GetNbins();
 305     double mK0SJac = 0.497614;            // mass of K0S (GeV)
 306     double nK0s = 1000;
 307     double deltay =0;
 308     double deltaeta = 0;
 309
 310
 311     TF1 *jacobian = new TF1("jacobian","[0]/TMath::Sqrt([1]*[1]*TMath::Power(TMath::CosH(x),2)+[0]*[0])",-etarange,etarange);
 312     jacobian->SetParameter(1,0.493);
 313     //here we're going to come up with a simple model that looks at a non-flat eta dependence.
 314     //if I assume that the pseudorapidity distribution has the shape y = m|x|+b, I have to normalize by the integral of that function after I take the jacobian so I can get the right weighting.  That integral works out to be
 315     float m = +0.15;
 316     float b = 1;
 317     float renormalization = (m*etarange*etarange)/2+b*etarange;
 318     TF1 *jacobianReweighted = new TF1("jacobianReweighted","([2]*x+[3])*[4]*[0]/TMath::Sqrt([1]*[1]*TMath::Power(TMath::CosH(x),2)+[0]*[0])",-etarange,etarange);
 319     jacobianReweighted->SetParameter(1,0.493);
 320     jacobianReweighted->SetParameter(2,m);
 321     jacobianReweighted->SetParameter(3,b);
 322     jacobianReweighted->SetParameter(4,renormalization);
 323     //==================CAIO'S TOY MODEL FOR JACOBIAN==========================
 324     float overallAverageJac[5] = {0,0,0,0,0};
 325     int totK0S = 0;
 326     for(int a = 1; a<NbinsJac; a ++){
 327
 328       int count = 0;
 329       double binContent = histoSys[centbin]->GetBinContent(a);
 330       double binError = histoSys[centbin]->GetBinError(a);
 331       jacobian->SetParameter(0,histoSys[centbin]->GetXaxis()->GetBinCenter(a));
 332       float averageJac = jacobian->Integral(-etarange,etarange) / (2*etarange);
 333       jacobian->SetParameter(0,histoSys[centbin]->GetXaxis()->GetBinLowEdge(a));
 334       float lowJac = jacobian->Integral(-etarange,etarange) / (2*etarange);
 335       jacobian->SetParameter(0,histoSys[centbin]->GetXaxis()->GetBinLowEdge(a+1));
 336       float highJac = jacobian->Integral(-etarange,etarange) / (2*etarange);
 337       float eta0Jac = jacobian->Eval(0.0);
 338       jacobianReweighted->SetParameter(0,histoSys[centbin]->GetXaxis()->GetBinCenter(a));
 339       float altJac = jacobianReweighted->Integral(-etarange,etarange);
 340       double correc =1;
 341
 342       double mult = (binContent*correc)*deltay*deltaeta;
 343       float scale = averageJac;
 344       switch(whichjacobian){
 345       case 1:
 346         scale = averageJac;
 347       case 2:
 348         scale = lowJac;
 349       case 3:
 350         scale = highJac;
 351       case 4:
 352         scale = eta0Jac;
 353       case 5:
 354         scale = altJac;
 355       }
 356       histoSpectrumK0SJac->SetBinContent(a,scale*binContent);
 357       histoSpectrumK0SJac->SetBinError(a,scale*binError);
 358       totK0S += binContent;
 359       overallAverageJac[0] += averageJac*binContent;
 360       overallAverageJac[1] += lowJac*binContent;
 361       overallAverageJac[2] += highJac*binContent;
 362       overallAverageJac[3] += eta0Jac*binContent;
 363       overallAverageJac[4] += altJac*binContent;
 364
 365     }
 366     cEmpty->cd();
 367     histoSys[centbin]->Draw();
 368     histoSpectrumK0SJac->Draw("same");
 369     cEmpty->SaveAs(Form("/tmp/KaonCut%i%s.png",etCutNum,det));
 370     float minJacobian = overallAverageJac[0]/totK0S;
 371     float maxJacobian = overallAverageJac[0]/totK0S;
 372     cout<<"Average jacobians: ";
 373     for(int i=0;i<=3;i++){
 374       cout<<" "<< overallAverageJac[i]/totK0S;
 375       if(minJacobian>overallAverageJac[i]/totK0S) minJacobian = overallAverageJac[i]/totK0S;
 376       if(maxJacobian<overallAverageJac[i]/totK0S) maxJacobian = overallAverageJac[i]/totK0S;
 377     }
 378     cout<<endl;
 379     float errJacobian = (maxJacobian-minJacobian)/2.0;
 380     float meanJacobian = (maxJacobian+minJacobian)/2.0;
 381     cout<<"Jacobian with error "<<meanJacobian<<" +/- "<<errJacobian<<endl;
 382     cEmpty->cd();
 383     funcBlastk0Jac[centbin] = FitParticle(histoSpectrumK0SJac,"K0",-1,-1,-1,0,0,3,-1,kFitBlastWave);
 384     funcBlastk0[centbin]->Draw("same");
 385
 386     // gets integral in input K0S histograms for different centrality bins
 387     fitINTEGRALdata[centbin] = histoSpectrumK0SJac->Integral("width");
 388
 389     cout<<"integral of fit for centrality bin "<<centbin<<": "<<fitINTEGRALdata[centbin]<<endl;
 390     cout<<funcBlastk0Jac[centbin]->Integral(0.0,100.0)<<endl;
 391
 392     //==================END CAIO'S TOY MODEL FOR JACOBIAN==========================
 393
 394     //==================THROWING RANDOM K0S==========================
 395     float totetcK0SOutOfAcc = 0;
 396     float totetcK0SOutOfAccCorr = 0;
 397     float avgCorrOutOfAcc = 0;
 398
 399     //THROW random K0S particle
 400     for(int i=0; i<nParticles; i++){
 401       float mK0S = 0.497614;            // mass of K0S (GeV)
 402       //You can change funBlastk0 to funBlastk0up or funBlastk0low for systematica error analysis
 403       float ptK0S = funcBlastk0Jac[centbin]->GetRandom(); // array? centbin =1
 404       histoSpectrumK0S->Fill(ptK0S);
 405
 406       // *******************************************
 407       float etK0S = sqrt(ptK0S*ptK0S + mK0S*mK0S);          // calculate transverse energy of K0S
 408       histoETspectrumK0S->Fill(etK0S);
 409       totET += etK0S;
 410
 411       int mybinK0S = my2DhistoK0S->GetXaxis()->FindBin(ptK0S);
 412       int xbins = my2DhistoK0S->GetNbinsX();
 413
 414       float scaleForOutOfAcc = hRatioOutOfAccOverInAcc->GetBinContent(hRatioOutOfAccOverInAcc->FindBin(ptK0S));
 415       //cout<<"Scale for out of acc "<<scaleForOutOfAcc<<endl;
 416
 417       float etcK0S = 0;
 418       float etcK0SOutOfAcc = 0;
 419       // get ET projection of K0S and sum
 420       if(mybinK0S>0 && mybinK0S<=bincount){  //nbins  // then this is in our range
 421         etcK0S = ((TH1D*)histoarrayK0S.At(mybinK0S-1))->GetRandom();
 422         etcK0SOutOfAcc = ((TH1D*)histoarrayK0SOutOfAcceptance.At(mybinK0S-1))->GetRandom();
 423         histoETCspectrumK0S->Fill(etcK0S);
 424         //here we deal with a problem.  The get random function gets confused for zeros.  It gives a small but non-zero value.  But we can never have a deposit less than the minimum energy!
 425         //this doesn't mess up the kaons in the
 426         if(etCutNum>0 && etK0S<etCutOffs[etCutNum-1]){etK0S = 0.0;}
 427         if(etCutNum>0 && etcK0SOutOfAcc<etCutOffs[etCutNum-1]){etcK0SOutOfAcc = 0.0;}
 428
 429         //if(etcK0SOutOfAcc>0){cout<<"et dep "<<etcK0SOutOfAcc<<" scale "<<scaleForOutOfAcc<<" pT "<<ptK0S<<endl;}
 430         //else{cout<<"No energy deposited"<<endl;}
 431         //else{cout<<"deposit is not small "<<etcK0SOutOfAcc<<endl;}
 432       }
 433       else{cerr<<"Did not find pt bin!"<<endl;}    // this should never be reached
 434       totetcK0S += etcK0S;
 435       totetcK0SOutOfAccCorr += etcK0SOutOfAcc*scaleForOutOfAcc;
 436       //if(etcK0SOutOfAcc>1e-3) cout<<"Deposit is not small "<<etcK0SOutOfAcc<<" corr "<<scaleForOutOfAcc<<" product "<<totetcK0SOutOfAccCorr<<" total "<<endl;
 437       totetcK0SOutOfAcc +=etcK0SOutOfAcc;
 438       avgCorrOutOfAcc += scaleForOutOfAcc;
 439     }
 440
 441     // print Calculated and Data projection values of scaled Et for K0S
 442     cout<<"Total ET observed (before renormalization): "<<(totetcK0S)<<endl;//Total ET observed
 443     cout<<"Total possible  energy (before renormalization):  "<<totET<<endl;
 444     cout<<"Average ET observed (simple): "<<(totetcK0S)/nParticles<<endl;
 445     //cout<<"Average ET observed out of acceptance (simple): "<<(totetcK0SOutOfAcc)/nParticles<<endl;
 446     //cout<<"Average corr for out of acc (simple): "<<(avgCorrOutOfAcc)/nParticles<<endl;
 447     //cout<<"ET observed out of acceptance corr (simple): "<<(totetcK0SOutOfAccCorr)<<endl;
 448     cout<<"Average ET observed out of acceptance corr (simple): "<<(totetcK0SOutOfAccCorr)/nParticles<<endl;
 449
 450     //Estimating the difference between the number of particles in |eta|<0.5 and our eta range
 451     //assume the distribution of particles is approximately N=m*|eta|+b
 452     //if N(eta=0) = 1, b=1
 453     //The percentage drop/increase in eta over one unit is m
 454     //The integral from 0-0.5 is
 455     m = 0.15;
 456     b = 1;
 457     float etaspectrameas = 0.5;
 458     float integralSpectraMeas = m*etaspectrameas*etaspectrameas/2+b*etaspectrameas;
 459     //and the integral over our range is:
 460     float integralMeasHigh = m*etarange*etarange/2+b*etarange;
 461     m = -0.15;
 462     float integralMeasLow = m*etarange*etarange/2+b*etarange;
 463     float etaRangeWeightMean = (integralMeasHigh/integralSpectraMeas+integralMeasLow/integralSpectraMeas)/2.0;
 464     float etaRangeWeightErr = (integralMeasHigh/integralSpectraMeas-integralMeasLow/integralSpectraMeas)/2.0;
 465     //cout<<" range high "<< integralMeasHigh/integralSpectraMeas<< " range low "<< integralMeasLow/integralSpectraMeas<<endl;
 466
 467     //cout<<"Eta range scale: "<<etaRangeWeightMean<<" +/- "<<etaRangeWeightErr<<endl;
 468
 469     //getting the final ET:
 470     //multiply by the dN/dy from the spectra paper (with error)
 471     //scale by Jacobian with the error on the Jacobian (with error)
 472     //multiply by the 4 kaon species
 473     //add an error for the fact that we need to extrapolate to our eta range
 474     float scale = 4*kaonYield[centbin]*meanJacobian*etaRangeWeightMean;
 475     float meanETperkaon = (totetcK0S+totetcK0SOutOfAccCorr)/nParticles;
 476     float meanETperkaonInAcc = (totetcK0S)/nParticles;
 477     float meanETperkaonOutOfAcc = (totetcK0SOutOfAccCorr)/nParticles;
 478     //now we're going to allow the ratio of kaons in and out of acceptance to vary by another 20%
 479     float fracerrFromInOutVariance = 0.2*meanETperkaonOutOfAcc/(meanETperkaonInAcc+meanETperkaonOutOfAcc);
 480     float fracerr = TMath::Sqrt(TMath::Power(kaonFracErr[centbin],2)+TMath::Power(errJacobian/meanJacobian,2)+TMath::Power(etaRangeWeightErr/etaRangeWeightMean,2));
 481     float fracerrTotal = TMath::Sqrt(TMath::Power(kaonFracErr[centbin],2)+TMath::Power(errJacobian/meanJacobian,2)+TMath::Power(etaRangeWeightErr/etaRangeWeightMean,2)+TMath::Power(fracerrFromInOutVariance,2));
 482     cout<<"Errors : Yield: "<<kaonFracErr[centbin]<<" Jacobian "<<errJacobian/meanJacobian<<" eta range "<<etaRangeWeightErr/etaRangeWeightMean<<" InOutVariance "<<fracerrFromInOutVariance<<endl;
 483     cout<<"Total ET deposited in acc: "<<meanETperkaonInAcc*scale<<" +/- "<<meanETperkaonInAcc*fracerr*scale<<endl;
 484     cout<<"Total ET deposited out of acc: "<<meanETperkaonOutOfAcc*scale<<" +/- "<<meanETperkaonOutOfAcc*fracerr*scale<<endl;
 485     cout<<"Total ET deposited: "<<meanETperkaon*scale<<" +/- "<<meanETperkaon*fracerrTotal*scale<<endl;
 486
 487     kaonDeposits[0][centbin] = meanETperkaon*scale;
 488     kaonDeposits[1][centbin] = meanETperkaon*scale*fracerrTotal;
 489
 490     delete histoSpectrumK0S;
 491     delete histoETspectrumK0S;
 492     delete histoETCspectrumK0S;
 493     delete histoSpectrumK0SJac;
 494     delete histoSpectrumK0SDivide;
 495     delete jacobian;
 496     delete jacobianReweighted;
 497   }
 498
 499
 500   PrintLatex(etCutNum,det,effCorr);
 501 }
 502
 503
 504 // function to define canvas characteristics globally
 505 void SetCanvasStyle(TCanvas *c){
 506   c->SetBorderSize(0);
 507   c->SetFillColor(0);
 508   c->SetBorderMode(0);
 509   c->SetFrameFillColor(0);
 510   c->SetFrameBorderMode(0);
 511   c->SetTopMargin(0.0254237);    // 0.04
 512   c->SetRightMargin(0.0322581);  // 0.04
 513   c->SetLeftMargin(0.129032);    // 0.181452
 514   c->SetBottomMargin(0.134409);
 515 }
 516
 517 // function to define characteristics globally of plots
 518 void SetStyles(TH1 *histo, char *ytitle, char *xtitle, int color, int marker){
 519   histo->Sumw2();
 520   histo->GetYaxis()->SetTitle(ytitle);
 521   histo->GetXaxis()->SetTitle(xtitle);
 522   histo->SetMarkerColor(color);
 523   histo->SetLineColor(color);
 524   histo->SetMarkerStyle(marker);
 525 }
 526
 527
 528 TLegend *GetLegend(float x1, float y1, float x2, float y2){
 529   TLegend *leg = new TLegend(x1,y1,x2,y2);
 530   leg->SetBorderSize(0);
 531   leg->SetTextFont(62);
 532   leg->SetLineColor(1);
 533   leg->SetLineStyle(1);
 534   leg->SetLineWidth(1);
 535   leg->SetFillColor(0);
 536   leg->SetFillStyle(0);
 537   leg->SetTextSize(0.0381356);
 538   return leg;
 539 }
 540
 541 void Draw(TH1 *data, TF1 *fLevy, TF1 *fBlast, int centbin,char *pid,char *name){
 542   TCanvas *canvas = new TCanvas("canvas","canvas",500,500);
 543   SetCanvasStyle(canvas);
 544   //canvas->SetLogy();
 545   data->SetTitle("p_{T} spectra of K$^{0}_{S}$ centbin=%i");
 546   data->GetYaxis()->SetTitle("1/2#pi p_{T} d^2N/dp_{T}dy");
 547   data->GetXaxis()->SetTitle("p_{T} (GeV)");
 548   data->GetXaxis()->SetTitleSize(0.05);
 549   data->GetYaxis()->SetTitleSize(0.05);
 550   data->GetXaxis()->SetRange(1,data->FindBin(3.5));
 551   //data->SetMarkerStyle(20);
 552   data->Draw();
 553   fLevy->Draw("same");
 554   fBlast->Draw("same");
 555   fLevy->SetLineColor(2);
 556   fBlast->SetLineColor(4);
 557   TLegend *legend = GetLegend(0.118952,0.120763,0.368952,0.271186);
 558   legend->AddEntry(data,Form("%s centrality bin %i",pid,centbin),"l");
 559   legend->AddEntry(fLevy,"Levy Fit");
 560   legend->AddEntry(fBlast,"BlastWave Fit");
 561   legend->Draw();
 562   canvas->SaveAs(Form("pics/%s%i.png",name,centbin));
 563   delete canvas;
 564 }
 565
 566 void PrintLatex(int etCutNum, char *det,Bool_t effCorr){
 567   TString cutNumString = Form("%i",etCutNum);
 568   ofstream myfile;
 569   TString tag = "";
 570   if(!effCorr){tag = "NoEffCorr";}
 571   TString tmpName = Form("datatablesKaonCut%i%s.tex",etCutNum,det);
 572   myfile.open(tmpName);
 573   cout<<"Creating "<<tmpName<<endl;
 574
 575   myfile<<Form("%3.2f",etCutOffs[etCutNum-1]);
 576   myfile<<"& ";
 577   //myfile<<"K$^{0}_{S}$  ";
 578   for(int j=0;j<nbins;j++){
 579     myfile<<Form("%3.1f $\\pm$ %3.1f",kaonDeposits[0][j],kaonDeposits[1][j]);
 580     if(j!=nbins-1) myfile<<"& ";
 581   }
 582   //myfile<<" & ";
 583   myfile<<"\\\\%data"<<endl;
 584   myfile<<endl<<endl;
 585
 586   myfile.close();
 587   cout<<"kaonCorr["<<nbins<<"] = {";
 588   for(int j=0;j<nbins;j++){
 589     cout<<kaonDeposits[0][j];
 590     if(j!=nbins-1) cout<<",";
 591   }
 592   cout<<"};"<<endl;
 593   cout<<"kaonError["<<nbins<<"] = {";
 594   for(int j=0;j<nbins;j++){
 595     cout<<kaonDeposits[1][j];
 596     if(j!=nbins-1) cout<<",";
 597   }
 598   cout<<"};"<<endl;
 599   cerr<<"I got here 587"<<endl;
 600   ofstream myfile2;
 601   cerr<<"I got here 589"<<endl;
 602   cerr<<"cut num "<<etCutNum<<" string "<<cutNumString<<endl;
 603   TString textfilename2 = Form("KaonCut%i%s%s.dat",etCutNum,det,tag.Data());
 604   //TString textfilename2 = "Kaons"+det+Form("%i",etCutNum)+".dat";
 605   cerr<<"I got here 590"<<endl;
 606   cout<<"Creating "<<textfilename2<<endl;
 607   myfile2.open (textfilename2.Data());
 608   for(int j=0;j<nbins;j++){
 609     myfile2<<Form("%2.3f %2.3f",kaonDeposits[0][j],kaonDeposits[1][j])<<endl;
 610   }
 611   myfile2.close();
 612 //   ofstream myfile3;
 613 //   TString textfilename2 = Form("KaonCut%i%sShort.dat",etCutNum,det);
 614 //   //TString textfilename = "Kaons"+det+Form("CutNum%i",etCutNum)+"Short.dat";
 615 //   cout<<"Creating "<<textfilename<<endl;
 616 //   myfile3.open (textfilename.Data());
 617 //   int startbin = 0;
 618 //   for(int j=0;j<10;j++){
 619 //     float mean = 0;
 620 //     float err = 0;
 621 //     if(j>1){
 622 //       startbin+=2;
 623 //       mean = (kaonDeposits[0][startbin]+kaonDeposits[0][startbin+1])/2.0;
 624 //       err = (kaonDeposits[1][startbin]+kaonDeposits[1][startbin+1])/2.0;
 625 //     }
 626 //     else{
 627 //       mean = kaonDeposits[0][startbin];
 628 //       err = kaonDeposits[1][startbin];
 629 //       startbin++;
 630 //     }
 631 //     myfile3<<Form("%2.3f %2.3f",kaonDeposits[0][j],kaonDeposits[1][j])<<endl;
 632 //   }
 633 //   myfile3.close();
 634
 635 }
 636
 637 void PrintArrays(){
 638   ofstream myfile;
 639   myfile.open("arraysv0.dat");
 640
 641   myfile<<"etk0[2][10] = {";
 642   for(int i=0;i<2;i++){
 643     myfile<<"{";
 644     for(int j=0;j<nbins;j++){
 645       myfile<<etk0[i][j];
 646       if(j<nbins-1) myfile<<",";
 647     }
 648     myfile<<"}";
 649     if(i<1) myfile<<",";
 650   }
 651   myfile<<"};"<<endl;
 652   myfile<<endl<<endl;
 653
 654   myfile<<"etBlastk0[2][10] = {";
 655   for(int i=0;i<2;i++){
 656     myfile<<"{";
 657     for(int j=0;j<nbins;j++){
 658       myfile<<etBlastk0[i][j];
 659       if(j<nbins-1) myfile<<",";
 660     }
 661     myfile<<"}";
 662     if(i<1) myfile<<",";
 663   }
 664   myfile<<"};"<<endl;
 665   myfile<<endl<<endl;
 666
 667   myfile<<"etBlastk02D[2][10] = {";
 668   for(int i=0;i<2;i++){
 669     myfile<<"{";
 670     for(int j=0;j<nbins;j++){
 671       myfile<<etBlastk02D[i][j];
 672       if(j<nbins-1) myfile<<",";
 673     }
 674     myfile<<"}";
 675     if(i<1) myfile<<",";
 676   }
 677   myfile<<"};"<<endl;
 678
 679   myfile.close();
 680 }
 681
 682
 683
 684
 685
 686
 687