.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / 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");
44   gSystem->Load("libVMC");
45   gSystem->Load("libMinuit");
46   gSystem->Load("libSTEERBase");
47   gSystem->Load("libESD");
48   gSystem->Load("libAOD");
49   gSystem->Load("libANALYSIS");
50   gSystem->Load("libANALYSISalice");
51   gSystem->Load("libCORRFW");
52   gSystem->Load("libPWGLFspectra");
53   gSystem->Load("libPWGTools");
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