#include <TObjArray.h>
#include <TGraphErrors.h>
#include <TString.h>
+#include <TSpline.h>
+#include <TRandom3.h>
#include "AliVParticle.h"
#include "AliMCParticle.h"
}//resonance cut
// HBT like cut
- if(fHBTCut){ // VERSION 3 (all pairs)
- //if(fHBTCut && charge1 * charge2 > 0){ // VERSION 2 (only for LS)
+ //if(fHBTCut){ // VERSION 3 (all pairs)
+ if(fHBTCut && charge1 * charge2 > 0){ // VERSION 2 (only for LS)
//if( dphi < 3 || deta < 0.01 ){ // VERSION 1
// continue;
Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
- const Float_t kLimit = 0.02 * 3;
+ const Float_t kLimit = fHBTCutValue * 3;
Float_t dphistarminabs = 1e5;
- Float_t dphistarmin = 1e5;
+ //Float_t dphistarmin = 1e5;
if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
Float_t dphistarabs = TMath::Abs(dphistar);
if (dphistarabs < dphistarminabs) {
- dphistarmin = dphistar;
+ //dphistarmin = dphistar;
dphistarminabs = dphistarabs;
}
}
- if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02) {
+ if (dphistarminabs < fHBTCutValue && TMath::Abs(deta) < fHBTCutValue) {
//AliInfo(Form("HBT: Removed track pair %d %d with [[%f %f]] %f %f %f | %f %f %d %f %f %d %f", i, j, deta, dphi, dphistarminabs, dphistar1, dphistar2, phi1rad, pt1, charge1, phi2rad, pt2, charge2, bSign));
continue;
}
hTemp3->Sumw2();
hTemp4->Sumw2();
hTemp1->Add(hTemp3,-1.);
- hTemp1->Scale(1./hTemp5->GetEntries());
+ hTemp1->Scale(1./hTemp5->Integral());
hTemp2->Add(hTemp4,-1.);
- hTemp2->Scale(1./hTemp6->GetEntries());
+ hTemp2->Scale(1./hTemp6->Integral());
gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
gHistBalanceFunctionHistogram->Scale(0.5);
//Printf("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0));
// Project into the wanted space (1st: analysis step, 2nd: axis)
- TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
- TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
- TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
- TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
+ TH1D* hTempHelper1 = (TH1D*)fHistPN->Project(0,iVariablePair);
+ TH1D* hTempHelper2 = (TH1D*)fHistNP->Project(0,iVariablePair);
+ TH1D* hTempHelper3 = (TH1D*)fHistPP->Project(0,iVariablePair);
+ TH1D* hTempHelper4 = (TH1D*)fHistNN->Project(0,iVariablePair);
TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
// ============================================================================================
// the same for event mixing
- TH1D* hTemp1Mix = (TH1D*)fHistPNMix->Project(0,iVariablePair);
- TH1D* hTemp2Mix = (TH1D*)fHistNPMix->Project(0,iVariablePair);
- TH1D* hTemp3Mix = (TH1D*)fHistPPMix->Project(0,iVariablePair);
- TH1D* hTemp4Mix = (TH1D*)fHistNNMix->Project(0,iVariablePair);
+ TH1D* hTempHelper1Mix = (TH1D*)fHistPNMix->Project(0,iVariablePair);
+ TH1D* hTempHelper2Mix = (TH1D*)fHistNPMix->Project(0,iVariablePair);
+ TH1D* hTempHelper3Mix = (TH1D*)fHistPPMix->Project(0,iVariablePair);
+ TH1D* hTempHelper4Mix = (TH1D*)fHistNNMix->Project(0,iVariablePair);
TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
// ============================================================================================
-
+
+ hTempHelper1->Sumw2();
+ hTempHelper2->Sumw2();
+ hTempHelper3->Sumw2();
+ hTempHelper4->Sumw2();
+ hTemp5->Sumw2();
+ hTemp6->Sumw2();
+ hTempHelper1Mix->Sumw2();
+ hTempHelper2Mix->Sumw2();
+ hTempHelper3Mix->Sumw2();
+ hTempHelper4Mix->Sumw2();
+ hTemp5Mix->Sumw2();
+ hTemp6Mix->Sumw2();
+
+ // first put everything on positive x - axis (to be comparable with published data) --> ONLY IN THIS OPTION!
+
+ Double_t helperEndBin = 1.6;
+ if(iVariablePair==2) helperEndBin = TMath::Pi();
+
+ TH1D* hTempPos1 = new TH1D(Form("hTempPos1_%d_%d",iBinPsi,iBinVertex),Form("hTempPos1_%d_%d",iBinPsi,iBinVertex),hTempHelper1->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTempPos2 = new TH1D(Form("hTempPos2_%d_%d",iBinPsi,iBinVertex),Form("hTempPos2_%d_%d",iBinPsi,iBinVertex),hTempHelper2->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTempPos3 = new TH1D(Form("hTempPos3_%d_%d",iBinPsi,iBinVertex),Form("hTempPos3_%d_%d",iBinPsi,iBinVertex),hTempHelper3->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTempPos4 = new TH1D(Form("hTempPos4_%d_%d",iBinPsi,iBinVertex),Form("hTempPos4_%d_%d",iBinPsi,iBinVertex),hTempHelper4->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTempPos1Mix = new TH1D(Form("hTempPos1Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos1Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper1Mix->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTempPos2Mix = new TH1D(Form("hTempPos2Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos2Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper2Mix->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTempPos3Mix = new TH1D(Form("hTempPos3Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos3Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper3Mix->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTempPos4Mix = new TH1D(Form("hTempPos4Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos4Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper4Mix->GetNbinsX()/2,0,helperEndBin);
+
+ TH1D* hTemp1 = new TH1D(Form("hTemp1_%d_%d",iBinPsi,iBinVertex),Form("hTemp1_%d_%d",iBinPsi,iBinVertex),hTempHelper1->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTemp2 = new TH1D(Form("hTemp2_%d_%d",iBinPsi,iBinVertex),Form("hTemp2_%d_%d",iBinPsi,iBinVertex),hTempHelper2->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTemp3 = new TH1D(Form("hTemp3_%d_%d",iBinPsi,iBinVertex),Form("hTemp3_%d_%d",iBinPsi,iBinVertex),hTempHelper3->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTemp4 = new TH1D(Form("hTemp4_%d_%d",iBinPsi,iBinVertex),Form("hTemp4_%d_%d",iBinPsi,iBinVertex),hTempHelper4->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTemp1Mix = new TH1D(Form("hTemp1Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp1Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper1Mix->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTemp2Mix = new TH1D(Form("hTemp2Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp2Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper2Mix->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTemp3Mix = new TH1D(Form("hTemp3Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp3Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper3Mix->GetNbinsX()/2,0,helperEndBin);
+ TH1D* hTemp4Mix = new TH1D(Form("hTemp4Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp4Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper4Mix->GetNbinsX()/2,0,helperEndBin);
+
+
+ hTempPos1->Sumw2();
+ hTempPos2->Sumw2();
+ hTempPos3->Sumw2();
+ hTempPos4->Sumw2();
+ hTempPos1Mix->Sumw2();
+ hTempPos2Mix->Sumw2();
+ hTempPos3Mix->Sumw2();
+ hTempPos4Mix->Sumw2();
+
hTemp1->Sumw2();
hTemp2->Sumw2();
hTemp3->Sumw2();
hTemp2Mix->Sumw2();
hTemp3Mix->Sumw2();
hTemp4Mix->Sumw2();
-
- hTemp1->Scale(1./hTemp5->GetEntries());
- hTemp3->Scale(1./hTemp5->GetEntries());
- hTemp2->Scale(1./hTemp6->GetEntries());
- hTemp4->Scale(1./hTemp6->GetEntries());
+
+
+ for(Int_t i=0;i<hTempHelper1->GetNbinsX();i++){
+
+ Double_t binCenter = hTempHelper1->GetXaxis()->GetBinCenter(i+1);
+
+ if(iVariablePair==1){
+ if(binCenter>0){
+ hTempPos1->SetBinContent(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinContent(i+1));
+ hTempPos2->SetBinContent(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinContent(i+1));
+ hTempPos3->SetBinContent(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinContent(i+1));
+ hTempPos4->SetBinContent(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinContent(i+1));
+ hTempPos1Mix->SetBinContent(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinContent(i+1));
+ hTempPos2Mix->SetBinContent(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinContent(i+1));
+ hTempPos3Mix->SetBinContent(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinContent(i+1));
+ hTempPos4Mix->SetBinContent(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinContent(i+1));
+
+ hTempPos1->SetBinError(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinError(i+1));
+ hTempPos2->SetBinError(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinError(i+1));
+ hTempPos3->SetBinError(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinError(i+1));
+ hTempPos4->SetBinError(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinError(i+1));
+ hTempPos1Mix->SetBinError(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinError(i+1));
+ hTempPos2Mix->SetBinError(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinError(i+1));
+ hTempPos3Mix->SetBinError(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinError(i+1));
+ hTempPos4Mix->SetBinError(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinError(i+1));
+ }
+ else{
+ hTemp1->SetBinContent(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinContent(i+1));
+ hTemp2->SetBinContent(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinContent(i+1));
+ hTemp3->SetBinContent(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinContent(i+1));
+ hTemp4->SetBinContent(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinContent(i+1));
+ hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinContent(i+1));
+ hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinContent(i+1));
+ hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinContent(i+1));
+ hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinContent(i+1));
+
+ hTemp1->SetBinError(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinError(i+1));
+ hTemp2->SetBinError(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinError(i+1));
+ hTemp3->SetBinError(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinError(i+1));
+ hTemp4->SetBinError(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinError(i+1));
+ hTemp1Mix->SetBinError(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinError(i+1));
+ hTemp2Mix->SetBinError(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinError(i+1));
+ hTemp3Mix->SetBinError(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinError(i+1));
+ hTemp4Mix->SetBinError(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinError(i+1));
+ }
+ }
+ else if(iVariablePair==2){
+ if(binCenter>0 && binCenter<TMath::Pi()){
+ hTempPos1->SetBinContent(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinContent(i+1));
+ hTempPos2->SetBinContent(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinContent(i+1));
+ hTempPos3->SetBinContent(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinContent(i+1));
+ hTempPos4->SetBinContent(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinContent(i+1));
+ hTempPos1Mix->SetBinContent(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinContent(i+1));
+ hTempPos2Mix->SetBinContent(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinContent(i+1));
+ hTempPos3Mix->SetBinContent(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinContent(i+1));
+ hTempPos4Mix->SetBinContent(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinContent(i+1));
+
+ hTempPos1->SetBinError(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinError(i+1));
+ hTempPos2->SetBinError(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinError(i+1));
+ hTempPos3->SetBinError(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinError(i+1));
+ hTempPos4->SetBinError(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinError(i+1));
+ hTempPos1Mix->SetBinError(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinError(i+1));
+ hTempPos2Mix->SetBinError(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinError(i+1));
+ hTempPos3Mix->SetBinError(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinError(i+1));
+ hTempPos4Mix->SetBinError(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinError(i+1));
+ }
+ else if(binCenter<0){
+ hTemp1->SetBinContent(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinContent(i+1));
+ hTemp2->SetBinContent(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinContent(i+1));
+ hTemp3->SetBinContent(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinContent(i+1));
+ hTemp4->SetBinContent(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinContent(i+1));
+ hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinContent(i+1));
+ hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinContent(i+1));
+ hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinContent(i+1));
+ hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinContent(i+1));
+
+ hTemp1->SetBinError(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinError(i+1));
+ hTemp2->SetBinError(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinError(i+1));
+ hTemp3->SetBinError(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinError(i+1));
+ hTemp4->SetBinError(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinError(i+1));
+ hTemp1Mix->SetBinError(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinError(i+1));
+ hTemp2Mix->SetBinError(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinError(i+1));
+ hTemp3Mix->SetBinError(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinError(i+1));
+ hTemp4Mix->SetBinError(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinError(i+1));
+ }
+ else{
+ hTemp1->SetBinContent(hTemp1->FindBin(TMath::Pi()-binCenter),hTempHelper1->GetBinContent(i+1));
+ hTemp2->SetBinContent(hTemp2->FindBin(TMath::Pi()-binCenter),hTempHelper2->GetBinContent(i+1));
+ hTemp3->SetBinContent(hTemp3->FindBin(TMath::Pi()-binCenter),hTempHelper3->GetBinContent(i+1));
+ hTemp4->SetBinContent(hTemp4->FindBin(TMath::Pi()-binCenter),hTempHelper4->GetBinContent(i+1));
+ hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(TMath::Pi()-binCenter),hTempHelper1Mix->GetBinContent(i+1));
+ hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(TMath::Pi()-binCenter),hTempHelper2Mix->GetBinContent(i+1));
+ hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(TMath::Pi()-binCenter),hTempHelper3Mix->GetBinContent(i+1));
+ hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(TMath::Pi()-binCenter),hTempHelper4Mix->GetBinContent(i+1));
+
+ hTemp1->SetBinError(hTemp1->FindBin(TMath::Pi()-binCenter),hTempHelper1->GetBinError(i+1));
+ hTemp2->SetBinError(hTemp2->FindBin(TMath::Pi()-binCenter),hTempHelper2->GetBinError(i+1));
+ hTemp3->SetBinError(hTemp3->FindBin(TMath::Pi()-binCenter),hTempHelper3->GetBinError(i+1));
+ hTemp4->SetBinError(hTemp4->FindBin(TMath::Pi()-binCenter),hTempHelper4->GetBinError(i+1));
+ hTemp1Mix->SetBinError(hTemp1Mix->FindBin(TMath::Pi()-binCenter),hTempHelper1Mix->GetBinError(i+1));
+ hTemp2Mix->SetBinError(hTemp2Mix->FindBin(TMath::Pi()-binCenter),hTempHelper2Mix->GetBinError(i+1));
+ hTemp3Mix->SetBinError(hTemp3Mix->FindBin(TMath::Pi()-binCenter),hTempHelper3Mix->GetBinError(i+1));
+ hTemp4Mix->SetBinError(hTemp4Mix->FindBin(TMath::Pi()-binCenter),hTempHelper4Mix->GetBinError(i+1));
+ }
+ }
+ }
+
+ hTemp1->Add(hTempPos1);
+ hTemp2->Add(hTempPos2);
+ hTemp3->Add(hTempPos3);
+ hTemp4->Add(hTempPos4);
+ hTemp1Mix->Add(hTempPos1Mix);
+ hTemp2Mix->Add(hTempPos2Mix);
+ hTemp3Mix->Add(hTempPos3Mix);
+ hTemp4Mix->Add(hTempPos4Mix);
+
+ hTemp1->Scale(1./hTemp5->Integral());
+ hTemp3->Scale(1./hTemp5->Integral());
+ hTemp2->Scale(1./hTemp6->Integral());
+ hTemp4->Scale(1./hTemp6->Integral());
// normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
// does not work here, so normalize also to trigger particles
// --> careful: gives different integrals then as with full 2D method
- hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
- hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
- hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
- hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
+ hTemp1Mix->Scale(1./hTemp5Mix->Integral());
+ hTemp3Mix->Scale(1./hTemp5Mix->Integral());
+ hTemp2Mix->Scale(1./hTemp6Mix->Integral());
+ hTemp4Mix->Scale(1./hTemp6Mix->Integral());
hTemp1->Divide(hTemp1Mix);
hTemp2->Divide(hTemp2Mix);
h4->Add(hTemp4);
}
+ delete hTemp1;
+ delete hTemp2;
+ delete hTemp3;
+ delete hTemp4;
+ delete hTemp1Mix;
+ delete hTemp2Mix;
+ delete hTemp3Mix;
+ delete hTemp4Mix;
+
+ delete hTempPos1;
+ delete hTempPos2;
+ delete hTempPos3;
+ delete hTempPos4;
+ delete hTempPos1Mix;
+ delete hTempPos2Mix;
+ delete hTempPos3Mix;
+ delete hTempPos4Mix;
}
}
hTemp3->Sumw2();
hTemp4->Sumw2();
hTemp1->Add(hTemp3,-1.);
- hTemp1->Scale(1./hTemp5->GetEntries());
+ hTemp1->Scale(1./hTemp5->Integral());
hTemp2->Add(hTemp4,-1.);
- hTemp2->Scale(1./hTemp6->GetEntries());
+ hTemp2->Scale(1./hTemp6->Integral());
gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
gHistBalanceFunctionHistogram->Scale(0.5);
hTemp3Mix->Sumw2();
hTemp4Mix->Sumw2();
- hTemp1->Scale(1./hTemp5->GetEntries());
- hTemp3->Scale(1./hTemp5->GetEntries());
- hTemp2->Scale(1./hTemp6->GetEntries());
- hTemp4->Scale(1./hTemp6->GetEntries());
+ hTemp1->Scale(1./hTemp5->Integral());
+ hTemp3->Scale(1./hTemp5->Integral());
+ hTemp2->Scale(1./hTemp6->Integral());
+ hTemp4->Scale(1./hTemp6->Integral());
// normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX());
binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
}
- Double_t binPsiLowEdge = 0.;
- Double_t binPsiUpEdge = 0.;
- Double_t binVertexLowEdge = 0.;
- Double_t binVertexUpEdge = 0.;
TH2D* h1 = NULL;
TH2D* h2 = NULL;
cout<<"In the balance function (2D->1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
- // determine the bin edges for this bin
- binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
- binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
- if(fVertexBinning){
- binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
- binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
- }
-
// Psi_2
fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
hTemp3Mix->Sumw2();
hTemp4Mix->Sumw2();
- hTemp1->Scale(1./hTemp5->GetEntries());
- hTemp3->Scale(1./hTemp5->GetEntries());
- hTemp2->Scale(1./hTemp6->GetEntries());
- hTemp4->Scale(1./hTemp6->GetEntries());
+ hTemp1->Scale(1./hTemp5->Integral());
+ hTemp3->Scale(1./hTemp5->Integral());
+ hTemp2->Scale(1./hTemp6->Integral());
+ hTemp4->Scale(1./hTemp6->Integral());
// normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX());
return gHistBalanceFunctionHistogram;
}
+//____________________________________________________________________//
+TH1D *AliBalancePsi::GetTriggers(TString type,
+ Double_t psiMin,
+ Double_t psiMax,
+ Double_t vertexZMin,
+ Double_t vertexZMax,
+ Double_t ptTriggerMin,
+ Double_t ptTriggerMax
+ ) {
+
+ // Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs,
+ // does the division by event mixing inside,
+ // and averaging over several vertexZ and centrality bins
+
+ // security checks
+ if(type != "PN" && type != "NP" && type != "PP" && type != "NN" && type != "ALL"){
+ AliError("Only types allowed: PN,NP,PP,NN,ALL");
+ return NULL;
+ }
+
+ TH1D *gHist = NULL;
+
+ // ranges (in bins) for vertexZ and centrality (psi)
+ Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin+0.00001);
+ Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
+ Int_t binVertexMin = 0;
+ Int_t binVertexMax = 0;
+ if(fVertexBinning){
+ binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin+0.00001);
+ binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
+ }
+ Double_t binPsiLowEdge = 0.;
+ Double_t binPsiUpEdge = 1000.;
+ Double_t binVertexLowEdge = 0.;
+ Double_t binVertexUpEdge = 1000.;
+
+
+ // first set to range and then obtain histogram
+ if(type=="PN" || type=="PP"){
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
+ gHist = (TH1D*)fHistP->Project(0,1);
+ }
+ else if(type=="NP" || type=="NN"){
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
+ gHist = (TH1D*)fHistN->Project(0,1);
+ }
+ else if(type=="ALL"){
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
+ gHist = (TH1D*)fHistN->Project(0,1);
+ gHist->Add((TH1D*)fHistP->Project(0,1));
+ }
+
+ return gHist;
+}
//____________________________________________________________________//
TH2D *AliBalancePsi::GetCorrelationFunction(TString type,
Double_t ptTriggerMax,
Double_t ptAssociatedMin,
Double_t ptAssociatedMax,
- AliBalancePsi *bMixed) {
+ AliBalancePsi *bMixed,
+ Bool_t normToTrig,
+ Double_t normalizationRangePhi) {
// Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs,
// does the division by event mixing inside,
TH2D *fMixed = NULL;
// ranges (in bins) for vertexZ and centrality (psi)
- Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
+ Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin+0.00001);
Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
Int_t binVertexMin = 0;
Int_t binVertexMax = 0;
if(fVertexBinning){
- binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
+ binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin+0.00001);
binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
}
Double_t binPsiLowEdge = 0.;
Double_t binVertexLowEdge = 0.;
Double_t binVertexUpEdge = 1000.;
+ // if event mixing is empty, we can not add that sub bin
+ // need to record the number of triggers for correct normalization
+ Double_t nTrigSubBinEmpty = 0.;
+
// loop over all bins
for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
- cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
+ //cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
// determine the bin edges for this bin
- binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
- binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
+ binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi) + 0.00001;
+ binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi) - 0.00001;
if(fVertexBinning){
- binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
- binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
+ binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex) + 0.00001;
+ binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex) - 0.00001;
}
// get the 2D histograms for the correct type
fSame = GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
fMixed = bMixed->GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
}
+
+ if(fMixed && normToTrig && fMixed->Integral()>0){
+
+ // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
+ // do it only on away-side (due to two-track cuts): pi +- pi/6.
+ Int_t binXmin = fMixed->GetXaxis()->FindBin(0-10e-5);
+ Int_t binXmax = fMixed->GetXaxis()->FindBin(0+10e-5);
+ Double_t binsX = (Double_t)(binXmax - binXmin + 1);
+ Int_t binYmin = fMixed->GetYaxis()->FindBin(TMath::Pi() - normalizationRangePhi);
+ Int_t binYmax = fMixed->GetYaxis()->FindBin(TMath::Pi() + normalizationRangePhi - 0.00001);
+ Double_t binsY = (Double_t)(binYmax - binYmin + 1);
+
+ Double_t mixedNorm = fMixed->Integral(binXmin,binXmax,binYmin,binYmax);
+ mixedNorm /= binsX * binsY;
+
+ // finite bin correction
+ Double_t binWidthEta = fMixed->GetXaxis()->GetBinWidth(fMixed->GetNbinsX());
+ Double_t maxEta = fMixed->GetXaxis()->GetBinUpEdge(fMixed->GetNbinsX());
+
+ Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
+ //Printf("Finite bin correction: %f", finiteBinCorrection);
+ mixedNorm /= finiteBinCorrection;
+
+ fMixed->Scale(1./mixedNorm);
+ }
if(fSame && fMixed){
// then get the correlation function (divide fSame/fmixed)
- fSame->Divide(fMixed);
-
- // NEW averaging:
+
+ if(fMixed->Integral()>0)
+ fSame->Divide(fMixed);
+
+ // averaging with number of triggers:
// average over number of triggers in each sub-bin
Double_t NTrigSubBin = 0;
if(type=="PN" || type=="PP")
- NTrigSubBin = (Double_t)(fHistP->Project(0,1)->GetEntries());
+ NTrigSubBin = (Double_t)(fHistP->Project(0,1)->Integral());
else if(type=="NP" || type=="NN")
- NTrigSubBin = (Double_t)(fHistN->Project(0,1)->GetEntries());
+ NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral());
+ else if(type=="ALL")
+ NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral());
fSame->Scale(NTrigSubBin);
- // for the first: clone
- if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){
- gHist = (TH2D*)fSame->Clone();
+ // only if event mixing has enough statistics
+ if(fMixed->Integral()>0){
+
+ // for the first: clone
+ if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){
+ gHist = (TH2D*)fSame->Clone();
+ }
+ else{ // otherwise: add for averaging
+ gHist->Add(fSame);
+ }
}
- else{ // otherwise: add for averaging
- gHist->Add(fSame);
+
+ // otherwise record the number of triggers for correct normalization
+ else{
+ nTrigSubBinEmpty += NTrigSubBin;
}
+
}
}
}
if(gHist){
- // OLD averaging:
- // average over number of bins nbinsVertex * nbinsPsi
- // gHist->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
-
- // NEW averaging:
- // average over number of triggers in each sub-bin
+ // averaging with number of triggers:
// first set to full range and then obtain number of all triggers
Double_t NTrigAll = 0;
if(type=="PN" || type=="PP"){
fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
- NTrigAll = (Double_t)(fHistP->Project(0,1)->GetEntries());
+ NTrigAll = (Double_t)(fHistP->Project(0,1)->Integral());
}
else if(type=="NP" || type=="NN"){
fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
- NTrigAll = (Double_t)(fHistN->Project(0,1)->GetEntries());
+ NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral());
+ }
+ else if(type=="ALL"){
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
+ fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
+ fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
+ NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral());
}
+
+ // subtract number of triggers with empty sub bins for correct normalization
+ NTrigAll -= nTrigSubBinEmpty;
+
gHist->Scale(1./NTrigAll);
-
}
return gHist;
//AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
//AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
//AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
+ //AliInfo(Form("Integral (1D): %lf",(Double_t)(fHistP->Project(0,1)->Integral())));
+ //AliInfo(Form("Integral (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->Integral())));
//TCanvas *c2 = new TCanvas("c2","");
//c2->cd();
//fHistPN->Project(0,1,2)->DrawCopy("colz");
- if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
- gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
+ if((Double_t)(fHistP->Project(0,1)->Integral())>0)
+ gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
//normalize to bin width
gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
//Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
//Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
- if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
- gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
+ if((Double_t)(fHistN->Project(0,1)->Integral())>0)
+ gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
//normalize to bin width
gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
//Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
//Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
- if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
- gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
+ if((Double_t)(fHistP->Project(0,1)->Integral())>0)
+ gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
//normalize to bin width
gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
//Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
//Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
- if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
- gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
+ if((Double_t)(fHistN->Project(0,1)->Integral())>0)
+ gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
//normalize to bin width
gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
}
// pt associated
- if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
+ if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)){
fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
-
+ }
+
//0:step, 1: Delta eta, 2: Delta phi
TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
if(!gHistNN){
gHistNN->Add(gHistPN);
// divide by sum of + and - triggers
- if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
- gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries() + fHistN->Project(0,1)->GetEntries()));
+ if((Double_t)(fHistN->Project(0,1)->Integral())>0 && (Double_t)(fHistP->Project(0,1)->Integral())>0)
+ gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral()));
//normalize to bin width
gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
}
//____________________________________________________________________//
-
Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist, Bool_t kUseZYAM,
Double_t &mean, Double_t &meanError,
Double_t &sigma, Double_t &sigmaError,
AliFatal(Form("Tag %s not found in %s", tag, configuration));
return 0;
}
+
+//____________________________________________________________________//
+Double_t AliBalancePsi::GetFWHM(Int_t gDeltaEtaPhi, TH1D* gHist,
+ Double_t &fwhm_spline, Double_t &fwhmError) {
+
+ // helper method to calculate the fwhm and errors of a TH1D
+ if(gHist){
+ Int_t repeat = 10000;
+
+ if (gDeltaEtaPhi == 1){
+ fwhm_spline = 0.;
+ fwhmError = 0.;
+ gHist->Smooth(4);
+
+ Double_t xmin = gHist->GetXaxis()->GetXmin();
+ Double_t xmax = gHist->GetXaxis()->GetXmax();
+
+ //+++++++FWHM TSpline+++++++++++++++++++++++++++//
+ TSpline3 *spline3 = new TSpline3(gHist,"b2e2",0,0);
+ spline3->SetLineColor(kGreen+2);
+ spline3->SetLineWidth(2);
+ //spline3->Draw("SAME");
+
+ Int_t Nbin = gHist->GetNbinsX();
+ Int_t bin_max_hist_y = gHist->GetMaximumBin();
+ Double_t bins_center_x = gHist->GetBinCenter(bin_max_hist_y);
+ Double_t max_spline = spline3->Eval(bins_center_x);
+ Double_t y_spline = -999;
+ Double_t y_spline_r = -999;
+ Double_t x_spline_l = -999;
+ Double_t x_spline_r = -999;
+
+ for (Int_t i= 0; i < bin_max_hist_y*1000; i++){
+ y_spline = spline3->Eval(xmin + i*(bins_center_x - xmin)/ (bin_max_hist_y*1000));
+ if (y_spline > max_spline/2){
+ x_spline_l = xmin + (i-1)*(bins_center_x - xmin)/(bin_max_hist_y*1000);
+ break;
+ }
+ }
+ for (Int_t j= 0; j < bin_max_hist_y*1000; j++){
+ y_spline_r = spline3->Eval(bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000));
+ if (y_spline_r < max_spline/2){
+ x_spline_r = bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000);
+ break;
+ }
+ }
+ fwhm_spline = x_spline_r - x_spline_l;
+ //+++++++FWHM TSpline++++++++++++++++++++++++//
+
+ //+++++++Error Calculation SPLINE++++++++++++//
+ Double_t dmeans,dsigmas;
+ TSpline3 *spline1;
+ Int_t bin_max_hist_y1;
+ Double_t bins_center_x1;
+ Double_t max_spline1;
+ Double_t y_spline1 = - 999.;
+ Double_t y_spline_r1 = - 999.;
+ Double_t x_spline_l1 = -999.;
+ Double_t x_spline_r1 = -999.;
+ Double_t fwhm_spline1 = 0.;
+ Double_t fwhm_T = 0.;
+
+ TH1F *hEmpty2 = new TH1F("hEmpty2","hEmpty2",Nbin,xmin,xmax);
+ TRandom3 *rgauss = new TRandom3(0);
+ TH1F *hists = new TH1F("hists","hists",1000,1.,3);
+
+ for (Int_t f=0; f<repeat; f++){ //100
+ for (Int_t h=0; h<Nbin+1; h++){
+ dmeans = gHist->GetBinContent(h);
+ dsigmas = gHist->GetBinError(h);
+ Double_t rgauss_point = rgauss->Gaus(dmeans,dsigmas);
+ hEmpty2->SetBinContent(h,rgauss_point);
+
+ //++++++++++++//
+ hEmpty2->Smooth(4);
+ //++++++++++++//
+ }
+ spline1 = new TSpline3(hEmpty2,"b2e2",0,0);
+ spline1->SetLineColor(kMagenta+1);
+ spline1->SetLineWidth(1);
+ //spline1->Draw("SAME");
+
+ bin_max_hist_y1 = hEmpty2->GetMaximumBin();
+ bins_center_x1 = hEmpty2->GetBinCenter(bin_max_hist_y1);
+ max_spline1 = spline1->Eval(bins_center_x1);
+
+ for (Int_t i= 0; i < bin_max_hist_y1*1000; i++){
+ y_spline1 = spline3->Eval(xmin + i*(bins_center_x1 - xmin)/ (bin_max_hist_y1*1000));
+ if (y_spline1 > max_spline1/2){
+ x_spline_l1 = xmin + (i-1)*(bins_center_x1 - xmin)/(bin_max_hist_y1*1000);
+ break;
+ }
+ }
+ for (Int_t j= 0; j < bin_max_hist_y1*1000; j++){
+ y_spline_r1 = spline3->Eval(bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000));
+ if (y_spline_r1 < max_spline1/2){
+ x_spline_r1 = bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000);
+ break;
+ }
+ }
+
+ fwhm_spline1 = x_spline_r1 - x_spline_l1;
+ hists->Fill(fwhm_spline1);
+ fwhm_T += (fwhm_spline - fwhm_spline1)*(fwhm_spline - fwhm_spline1);
+ }
+
+ fwhmError = TMath::Sqrt(TMath::Abs(fwhm_T/(repeat-1)));
+ //+++++++Error Calculation SPLINE+++++++++++//
+ }
+ else{
+ fwhm_spline = 0.;
+ fwhmError = 0.;
+ gHist->Smooth(4);
+
+ Double_t xmin = gHist->GetXaxis()->GetXmin();
+ Double_t xmax = gHist->GetXaxis()->GetXmax();
+
+ //+++++++FWHM TSpline+++++++++++++++++++++++++++//
+ TSpline3 *spline3 = new TSpline3(gHist,"b2e2",0,0);
+ spline3->SetLineColor(kGreen+2);
+ spline3->SetLineWidth(2);
+ //spline3->Draw("SAME");
+
+ Int_t Nbin = gHist->GetNbinsX();
+ Int_t bin_max_hist_y = gHist->GetMaximumBin();
+ Double_t bins_center_x = gHist->GetBinCenter(bin_max_hist_y);
+ Double_t max_spline = spline3->Eval(bins_center_x);
+
+ Int_t bin_min_hist_y = gHist->GetMinimumBin();
+ Double_t bins_center_xmin = gHist->GetBinCenter(bin_min_hist_y);
+ Double_t min_spline = spline3->Eval(bins_center_xmin);
+
+ Double_t y_spline = -999.;
+ Double_t y_spline_r = -999.;
+ Double_t x_spline_l = -999.;
+ Double_t x_spline_r = -999.;
+
+ for (Int_t i= 0; i < bin_max_hist_y*1000; i++){
+ y_spline = spline3->Eval(xmin + i*(bins_center_x - xmin)/ (bin_max_hist_y*1000));
+ if (y_spline > (max_spline+min_spline)/2){
+ x_spline_l = xmin + (i-1)*(bins_center_x - xmin)/(bin_max_hist_y*1000);
+ break;
+ }
+ }
+ for (Int_t j= 0; j < bin_max_hist_y*1000; j++){
+ y_spline_r = spline3->Eval(bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000));
+ if (y_spline_r < (max_spline+min_spline)/2){
+ x_spline_r = bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000);
+ break;
+ }
+ }
+ fwhm_spline = x_spline_r - x_spline_l;
+ //+++++++FWHM TSpline++++++++++++++++++++++++//
+
+ //+++++++Error Calculation SPLINE++++++++++++//
+ Double_t dmeans,dsigmas;
+ TSpline3 *spline1;
+ Int_t bin_max_hist_y1;
+ Double_t bins_center_x1;
+ Double_t max_spline1;
+ Double_t y_spline1 = -999.;
+ Double_t y_spline_r1 = -999.;
+ Double_t x_spline_l1 = -999.;
+ Double_t x_spline_r1 = -999.;
+ Double_t fwhm_spline1 = 0.;
+ Double_t fwhm_T = 0.;
+
+ TH1F *hEmpty2 = new TH1F("hEmpty2","hEmpty2",Nbin,xmin,xmax);
+ TRandom3 *rgauss = new TRandom3(0);
+ Int_t bin_min_hist_y1;
+ Double_t bins_center_xmin1,min_spline1;
+
+ for (Int_t f=0; f<repeat; f++){ //100
+ for (Int_t h=0; h<Nbin+1; h++){
+ dmeans = gHist->GetBinContent(h);
+ dsigmas = gHist->GetBinError(h);
+ Double_t rgauss_point = rgauss->Gaus(dmeans,dsigmas);
+ hEmpty2->SetBinContent(h,rgauss_point);
+
+ //++++++++++++//
+ hEmpty2->Smooth(4);
+ //++++++++++++//
+ }
+ spline1 = new TSpline3(hEmpty2,"b2e2",0,0);
+ spline1->SetLineColor(kMagenta+1);
+ spline1->SetLineWidth(1);
+ //spline1->Draw("SAME");
+
+ bin_max_hist_y1 = hEmpty2->GetMaximumBin();
+ bins_center_x1 = hEmpty2->GetBinCenter(bin_max_hist_y1);
+ max_spline1 = spline1->Eval(bins_center_x1);
+
+ bin_min_hist_y1 = hEmpty2->GetMinimumBin();
+ bins_center_xmin1 = hEmpty2->GetBinCenter(bin_min_hist_y1);
+ min_spline1 = spline1->Eval(bins_center_xmin1);
+
+ for (Int_t i= 0; i < bin_max_hist_y1*1000; i++){
+ y_spline1 = spline3->Eval(xmin + i*(bins_center_x1 - xmin)/ (bin_max_hist_y1*1000));
+ if (y_spline1 > (max_spline1+min_spline1)/2){
+ x_spline_l1 = xmin + (i-1)*(bins_center_x1 - xmin)/(bin_max_hist_y1*1000);
+ break;
+ }
+ }
+ for (Int_t j= 0; j < bin_max_hist_y1*1000; j++){
+ y_spline_r1 = spline3->Eval(bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000));
+ if (y_spline_r1 < (max_spline1+min_spline1)/2){
+ x_spline_r1 = bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000);
+ break;
+ }
+ }
+
+ fwhm_spline1 = x_spline_r1 - x_spline_l1;
+ fwhm_T += (fwhm_spline - fwhm_spline1)*(fwhm_spline - fwhm_spline1);
+ }
+ fwhmError = TMath::Sqrt(TMath::Abs(fwhm_T/(repeat-1)));
+ }
+ }
+ return fwhm_spline;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff