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);
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
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};
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}};
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);
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;
69 //centrality bin # (want to loop over 5 bins)
70 int centbin = 0; // can change later
71 TString detector = det;
73 if(detector.Contains("EMC")){
78 //==================READ IN DATA==========================
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
84 TCanvas *cChK = new TCanvas("cChK","ChargedKaons",700,500);
86 TH1D *histoskch[] = {NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL, NULL,NULL};
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");
92 ((TH1D*)histoskch[i])->Draw();
95 ((TH1D*)histoskch[i])->Draw("same");
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");
105 ((TH1D*)histoskch[i])->Draw();
108 ((TH1D*)histoskch[i])->Draw("same");
113 TCanvas *cEmpty = new TCanvas("cEmpty","Empty",700,500);
115 TString particleName = "K0";
116 TString particleLatex = "K^{0}_S";
117 if(kaonSelection==1){
119 particleLatex = "K^{+}";
121 if(kaonSelection==-1){
123 particleLatex = "K^{-}";
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};
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};
139 TH1 *histoSys[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL}; //K0S Clone histogram
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};
151 Double_t Abaryon = -1;
152 Double_t Aantibaryon = 1;
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};
163 TRandom* randy = new TRandom();
165 // integral of fit to data
166 float fitINTEGRALdata[10] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL};
168 if(kaonSelection!=0){//charged kaons
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));
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));
185 float NbinsSys = histoSysUp[i]->GetXaxis()->GetNbins();
187 for(int x = 1; x<=NbinsSys; x++){
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);
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);
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());
206 if(!histoskch[i]) cerr<<"Warning! Histogram does not exist!"<<endl;
207 funck0[i] = FitParticle((TH1D*)histoskch[i],particleName.Data());
209 Draw(histoskch[i],funck0[i],funcBlastk0[i],i,particleLatex.Data(),particleName.Data());
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
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;
224 TFile *f = TFile::Open(inputfilename, "READ");
225 TList *l = (TList*)f->Get("out1");
226 TH3F *fHistK0EDepositsVsPtInAcceptance;
227 TH3F *fHistK0EDepositsVsPtOutOfAcceptance;
229 fHistK0EDepositsVsPtInAcceptance =(TH3F*) l->FindObject("fHistK0EDepositsVsPtInAcceptance");
230 fHistK0EDepositsVsPtOutOfAcceptance =(TH3F*) l->FindObject("fHistK0EDepositsVsPtOutOfAcceptance");
233 fHistK0EDepositsVsPtInAcceptance =(TH3F*) l->FindObject("fHistK0EGammaVsPtInAcceptance");
234 fHistK0EDepositsVsPtOutOfAcceptance =(TH3F*) l->FindObject("fHistK0EGammaVsPtOutOfAcceptance");
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);
247 //find number of bins of Caio's histogram
249 bincount = my2DhistoK0S->GetXaxis()->GetNbins();
250 cout<<"Number of bins in histo: "<<bincount<<endl;
252 // create array to hold histograms
253 TClonesArray histoarrayK0S("TH1D",bincount); //(nbins);
254 TClonesArray histoarrayK0SOutOfAcceptance("TH1D",bincount); //(nbins);
256 //==================PROJECTIONS FOR INPUT ET DEPOSITS==========================
257 //loop over bins in x to define projections (nbins)
258 for(int i=0; i<bincount; i++){
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);
266 for(int j=0;j<nbins;j++){
268 cout<<"Working on centrality bin "<<centbin<<endl;
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);
280 float totET = 0; // total ET of K0S's calculated
281 float totetcK0S = 0; // total ET of K0S's from data
286 // declare edges of first and last bin
287 int mybin100 = my2DhistoK0S->GetXaxis()->FindBin(0.1001);
288 int mybin0 = my2DhistoK0S->GetXaxis()->FindBin(0.0001);
290 //Jacobian correction
292 TH1F *histoSpectrumK0SJac;
293 histoSpectrumK0SJac = (TH1F*)histoSys[centbin]->Clone("Jac");
294 SetStyles(histoSpectrumK0SJac,"dN/dp_{T}","p_{T} (GeV)", color1, marker1);
296 TH1F *histoSpectrumK0SDivide = (TH1F*)histoSys[centbin]->Clone("JacDiv");
299 // TH1F *histoeta = new TH1F("histoeta","",100,-1,1);
300 // SetStyles(histoeta,"N","#eta", color1, marker1);
302 //Jacobian Correction
304 float NbinsJac = histoSys[centbin]->GetXaxis()->GetNbins();
305 double mK0SJac = 0.497614; // mass of K0S (GeV)
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
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};
326 for(int a = 1; a<NbinsJac; a ++){
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);
342 double mult = (binContent*correc)*deltay*deltaeta;
343 float scale = averageJac;
344 switch(whichjacobian){
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;
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;
379 float errJacobian = (maxJacobian-minJacobian)/2.0;
380 float meanJacobian = (maxJacobian+minJacobian)/2.0;
381 cout<<"Jacobian with error "<<meanJacobian<<" +/- "<<errJacobian<<endl;
383 funcBlastk0Jac[centbin] = FitParticle(histoSpectrumK0SJac,"K0",-1,-1,-1,0,0,3,-1,kFitBlastWave);
384 funcBlastk0[centbin]->Draw("same");
386 // gets integral in input K0S histograms for different centrality bins
387 fitINTEGRALdata[centbin] = histoSpectrumK0SJac->Integral("width");
389 cout<<"integral of fit for centrality bin "<<centbin<<": "<<fitINTEGRALdata[centbin]<<endl;
390 cout<<funcBlastk0Jac[centbin]->Integral(0.0,100.0)<<endl;
392 //==================END CAIO'S TOY MODEL FOR JACOBIAN==========================
394 //==================THROWING RANDOM K0S==========================
395 float totetcK0SOutOfAcc = 0;
396 float totetcK0SOutOfAccCorr = 0;
397 float avgCorrOutOfAcc = 0;
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);
406 // *******************************************
407 float etK0S = sqrt(ptK0S*ptK0S + mK0S*mK0S); // calculate transverse energy of K0S
408 histoETspectrumK0S->Fill(etK0S);
411 int mybinK0S = my2DhistoK0S->GetXaxis()->FindBin(ptK0S);
412 int xbins = my2DhistoK0S->GetNbinsX();
414 float scaleForOutOfAcc = hRatioOutOfAccOverInAcc->GetBinContent(hRatioOutOfAccOverInAcc->FindBin(ptK0S));
415 //cout<<"Scale for out of acc "<<scaleForOutOfAcc<<endl;
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;}
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;}
433 else{cerr<<"Did not find pt bin!"<<endl;} // this should never be reached
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;
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;
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
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;
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;
467 //cout<<"Eta range scale: "<<etaRangeWeightMean<<" +/- "<<etaRangeWeightErr<<endl;
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;
487 kaonDeposits[0][centbin] = meanETperkaon*scale;
488 kaonDeposits[1][centbin] = meanETperkaon*scale*fracerrTotal;
490 delete histoSpectrumK0S;
491 delete histoETspectrumK0S;
492 delete histoETCspectrumK0S;
493 delete histoSpectrumK0SJac;
494 delete histoSpectrumK0SDivide;
496 delete jacobianReweighted;
500 PrintLatex(etCutNum,det,effCorr);
504 // function to define canvas characteristics globally
505 void SetCanvasStyle(TCanvas *c){
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);
517 // function to define characteristics globally of plots
518 void SetStyles(TH1 *histo, char *ytitle, char *xtitle, int color, int marker){
520 histo->GetYaxis()->SetTitle(ytitle);
521 histo->GetXaxis()->SetTitle(xtitle);
522 histo->SetMarkerColor(color);
523 histo->SetLineColor(color);
524 histo->SetMarkerStyle(marker);
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);
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);
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);
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");
562 canvas->SaveAs(Form("pics/%s%i.png",name,centbin));
566 void PrintLatex(int etCutNum, char *det,Bool_t effCorr){
567 TString cutNumString = Form("%i",etCutNum);
570 if(!effCorr){tag = "NoEffCorr";}
571 TString tmpName = Form("datatablesKaonCut%i%s.tex",etCutNum,det);
572 myfile.open(tmpName);
573 cout<<"Creating "<<tmpName<<endl;
575 myfile<<Form("%3.2f",etCutOffs[etCutNum-1]);
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<<"& ";
583 myfile<<"\\\\%data"<<endl;
587 cout<<"kaonCorr["<<nbins<<"] = {";
588 for(int j=0;j<nbins;j++){
589 cout<<kaonDeposits[0][j];
590 if(j!=nbins-1) cout<<",";
593 cout<<"kaonError["<<nbins<<"] = {";
594 for(int j=0;j<nbins;j++){
595 cout<<kaonDeposits[1][j];
596 if(j!=nbins-1) cout<<",";
599 cerr<<"I got here 587"<<endl;
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;
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());
618 // for(int j=0;j<10;j++){
623 // mean = (kaonDeposits[0][startbin]+kaonDeposits[0][startbin+1])/2.0;
624 // err = (kaonDeposits[1][startbin]+kaonDeposits[1][startbin+1])/2.0;
627 // mean = kaonDeposits[0][startbin];
628 // err = kaonDeposits[1][startbin];
631 // myfile3<<Form("%2.3f %2.3f",kaonDeposits[0][j],kaonDeposits[1][j])<<endl;
639 myfile.open("arraysv0.dat");
641 myfile<<"etk0[2][10] = {";
642 for(int i=0;i<2;i++){
644 for(int j=0;j<nbins;j++){
646 if(j<nbins-1) myfile<<",";
654 myfile<<"etBlastk0[2][10] = {";
655 for(int i=0;i<2;i++){
657 for(int j=0;j<nbins;j++){
658 myfile<<etBlastk0[i][j];
659 if(j<nbins-1) myfile<<",";
667 myfile<<"etBlastk02D[2][10] = {";
668 for(int i=0;i<2;i++){
670 for(int j=0;j<nbins;j++){
671 myfile<<etBlastk02D[i][j];
672 if(j<nbins-1) myfile<<",";