+
/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
//____________________________________________________________
AliHFEspectrum::AliHFEspectrum(const char *name):
TNamed(name, ""),
- fCFContainers(NULL),
+ fCFContainers(new TObjArray(kDataContainerV0+1)),
fTemporaryObjects(NULL),
fCorrelation(NULL),
fBackground(NULL),
fIPanaCharmBgSubtract(kFALSE),
fIPanaConversionBgSubtract(kFALSE),
fIPanaNonHFEBgSubtract(kFALSE),
- fNonHFEbgMethod2(kFALSE),
- fNonHFEmode(0),
+ fIPParameterizedEff(kFALSE),
+ fNonHFEsyst(0),
+ fBeauty2ndMethod(kFALSE),
+ fIPEffCombinedSamples(kTRUE),
fNbDimensions(1),
fNMCEvents(0),
- fNMCbgEvents(0),
fStepMC(-1),
fStepTrue(-1),
fStepData(-1),
fStepAfterCutsV0(-1),
fStepGuessedUnfolding(-1),
fNumberOfIterations(5),
- fChargeChoosen(-1),
+ fChargeChoosen(kAllCharge),
fNCentralityBinAtTheEnd(0),
+ fTestCentralityLow(-1),
+ fTestCentralityHigh(-1),
+ fFillMoreCorrelationMatrix(kFALSE),
fHadronEffbyIPcut(NULL),
- fConversionEff(NULL),
- fNonHFEEff(NULL),
+ fConversionEffbgc(NULL),
+ fNonHFEEffbgc(NULL),
+ fBSpectrum2ndMethod(NULL),
+ fkBeauty2ndMethodfilename(""),
fBeamType(0),
- fDebugLevel(0)
+ fEtaSyst(kTRUE),
+ fDebugLevel(0),
+ fWriteToFile(kFALSE),
+ fUnfoldBG(kFALSE)
{
//
// Default constructor
//
for(Int_t k = 0; k < 20; k++){
- fNEvents[k] = 0;
- fLowBoundaryCentralityBinAtTheEnd[k] = 0;
- fHighBoundaryCentralityBinAtTheEnd[k] = 0;
+ fNEvents[k] = 0;
+ fNMCbgEvents[k] = 0;
+ fLowBoundaryCentralityBinAtTheEnd[k] = 0;
+ fHighBoundaryCentralityBinAtTheEnd[k] = 0;
+ if(k<kCentrality)
+ {
+ fEfficiencyTOFPIDD[k] = 0;
+ fEfficiencyesdTOFPIDD[k] = 0;
+ fEfficiencyIPCharmD[k] = 0;
+ fEfficiencyIPBeautyD[k] = 0;
+ fEfficiencyIPBeautyesdD[k] = 0;
+ fEfficiencyIPConversionD[k] = 0;
+ fEfficiencyIPNonhfeD[k] = 0;
+
+ fConversionEff[k] = 0;
+ fNonHFEEff[k] = 0;
+ fCharmEff[k] = 0;
+ fBeautyEff[k] = 0;
+ fEfficiencyCharmSigD[k] = 0;
+ fEfficiencyBeautySigD[k] = 0;
+ fEfficiencyBeautySigesdD[k] = 0;
+ }
}
memset(fEtaRange, 0, sizeof(Double_t) * 2);
+ memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
+ memset(fConvSourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
+ memset(fNonHFESourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
}
//____________________________________________________________
// If no inclusive spectrum, then take only efficiency map for beauty electron
// and the appropriate correlation matrix
//
- Int_t kNdim = 3;
- if(fBeamType==0) kNdim=3;
- if(fBeamType==1) kNdim=4;
-
- Int_t dims[kNdim];
- // Get the requested format
- if(fBeamType==0)
- {
- // pp analysis
- switch(fNbDimensions){
- case 1: dims[0] = 0;
- break;
- case 2: for(Int_t i = 0; i < 2; i++) dims[i] = i;
- break;
- case 3: for(Int_t i = 0; i < 3; i++) dims[i] = i;
- break;
- default:
- AliError("Container with this number of dimensions not foreseen (yet)");
- return kFALSE;
- };
- }
-
- if(fBeamType==1)
- {
- // PbPb analysis; centrality as first dimension
- Int_t nbDimensions = fNbDimensions;
- fNbDimensions = fNbDimensions + 1;
- switch(nbDimensions){
- case 1:
- dims[0] = 5;
- dims[1] = 0;
- break;
- case 2:
- dims[0] = 5;
- for(Int_t i = 0; i < 2; i++) dims[(i+1)] = i;
- break;
- case 3:
- dims[0] = 5;
- for(Int_t i = 0; i < 3; i++) dims[(i+1)] = i;
- break;
- default:
- AliError("Container with this number of dimensions not foreseen (yet)");
- return kFALSE;
- };
- }
-
+
+ if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
+
+ Int_t kNdim = 3;
+ Int_t kNcentr =1;
+ Int_t ptpr =0;
+ if(fBeamType==0) kNdim=3;
+ if(fBeamType==1)
+ {
+ kNdim=4;
+ kNcentr=11;
+ ptpr=1;
+ }
+
+ Int_t dims[kNdim];
+ // Get the requested format
+ if(fBeamType==0){
+ // pp analysis
+ switch(fNbDimensions){
+ case 1: dims[0] = 0;
+ break;
+ case 2: for(Int_t i = 0; i < 2; i++) dims[i] = i;
+ break;
+ case 3: for(Int_t i = 0; i < 3; i++) dims[i] = i;
+ break;
+ default:
+ AliError("Container with this number of dimensions not foreseen (yet)");
+ return kFALSE;
+ };
+ }
+
+ if(fBeamType==1){
+ // PbPb analysis; centrality as first dimension
+ Int_t nbDimensions = fNbDimensions;
+ fNbDimensions = fNbDimensions + 1;
+ switch(nbDimensions){
+ case 1: dims[0] = 5;
+ dims[1] = 0;
+ break;
+ case 2: dims[0] = 5;
+ for(Int_t i = 0; i < 2; i++) dims[(i+1)] = i;
+ break;
+ case 3: dims[0] = 5;
+ for(Int_t i = 0; i < 3; i++) dims[(i+1)] = i;
+ break;
+ default:
+ AliError("Container with this number of dimensions not foreseen (yet)");
+ return kFALSE;
+ };
+ }
// Data container: raw spectrum + hadron contamination
AliCFContainer *datacontainer = 0x0;
datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
}
else{
- datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
+
+ datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
}
AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
- AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen);
- AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen);
+ AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+ AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
SetContainer(datacontainerD,AliHFEspectrum::kDataContainer);
// MC container: ESD/MC efficiency maps + MC/MC efficiency maps
AliCFContainer *mccontaineresd = 0x0;
+ AliCFContainer *mccontaineresdbg = 0x0;
AliCFContainer *mccontainermc = 0x0;
AliCFContainer *mccontainermcbg = 0x0;
AliCFContainer *nonHFEweightedContainer = 0x0;
AliCFContainer *convweightedContainer = 0x0;
+ AliCFContainer *nonHFEtempContainer = 0x0;//temporary container to be sliced for the fnonHFESourceContainers
+ AliCFContainer *convtempContainer = 0x0;//temporary container to be sliced for the fConvSourceContainers
if(fInclusiveSpectrum) {
mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
}
else {
mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
+ mccontaineresdbg = bghfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
- nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
- convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
- if((!nonHFEweightedContainer) || (!convweightedContainer)) return kFALSE;
+
+ if(fNonHFEsyst){
+ const Char_t *sourceName[kElecBgSources]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
+ const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
+ for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
+ for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
+ nonHFEtempContainer = bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]));
+ convtempContainer = bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
+ for(Int_t icentr=0;icentr<kNcentr;icentr++)
+ {
+ if(fBeamType==0)
+ {
+ fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
+ fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
+ }
+ if(fBeamType==1)
+ {
+
+ fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
+ fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
+ }
+// if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
+ }
+ if(fBeamType == 1)break;
+ }
+ }
+ }
+ // else{
+ nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
+ convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
+ if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE;
+ //}
}
if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
-
+
Int_t source = -1;
if(!fInclusiveSpectrum) source = 1; //beauty
- AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen);
- AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen);
+ AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+ AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC);
SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD);
mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC);
- source = 2; //conversion
- mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
- AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
- efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. be careful with step #
- fConversionEff = (TH1D*)efficiencyConv->Project(0);
+ //if(!fNonHFEsyst){
+ AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
+ SetContainer(nonHFEweightedContainerD,AliHFEspectrum::kMCWeightedContainerNonHFEESD);
+ AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
+ SetContainer(convweightedContainerD,AliHFEspectrum::kMCWeightedContainerConversionESD);
+ //}
- source = 3; //non HFE except for the conversion
- mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
- AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
- efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. be careful with step #
- fNonHFEEff = (TH1D*)efficiencyNonhfe->Project(0);
+ SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims);
- AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
- SetContainer(nonHFEweightedContainerD,AliHFEspectrum::kMCWeightedContainerNonHFEESD);
- AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
- SetContainer(convweightedContainerD,AliHFEspectrum::kMCWeightedContainerConversionESD);
}
-
-// MC container: correlation matrix
+ // MC container: correlation matrix
THnSparseF *mccorrelation = 0x0;
if(fInclusiveSpectrum) {
- if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
- else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
+ if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
+ else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
//else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
if(!mccorrelation) return kFALSE;
- THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims);
+ Int_t testCentralityLow = fTestCentralityLow;
+ Int_t testCentralityHigh = fTestCentralityHigh;
+ if(fFillMoreCorrelationMatrix) {
+ testCentralityLow = fTestCentralityLow-1;
+ testCentralityHigh = fTestCentralityHigh+1;
+ }
+ THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims,testCentralityLow,testCentralityHigh);
if(!mccorrelationD) {
printf("No correlation\n");
return kFALSE;
if(v0hfecontainer) {
AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
if(!containerV0) return kFALSE;
- AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron+1);
+ AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
if(!containerV0Electron) return kFALSE;
SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0);
}
ccorrelation->cd(1);
if(mccorrelationD) {
if(fBeamType==0){
- TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions,0);
- ptcorrelation->Draw("colz");
+ TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions,0);
+ ptcorrelation->Draw("colz");
}
if(fBeamType==1){
- TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions+1,1);
- ptcorrelation->Draw("colz");
+ TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions+1,1);
+ ptcorrelation->Draw("colz");
}
}
+ if(fWriteToFile){
+ ccontaminationspectrum->SaveAs("contaminationspectrum.eps");
+ ccorrelation->SaveAs("correlationMatrix.eps");
+ }
}
+ TFile *file = TFile::Open("tests.root","recreate");
+ datacontainerD->Write("data");
+ mccontainermcD->Write("mcefficiency");
+ mccorrelationD->Write("correlationmatrix");
+ file->Close();
return kTRUE;
-
-
+}
+
+
+//____________________________________________________________
+void AliHFEspectrum::CallInputFileForBeauty2ndMethod(){
+ //
+ // get spectrum for beauty 2nd method
+ //
+ //
+ TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename);
+ fBSpectrum2ndMethod = new TH1D(*static_cast<TH1D*>(inputfile2ndmethod->Get("BSpectrum")));
}
ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
ratiocorrected->SetStats(0);
ratiocorrected->Draw();
+ if(fWriteToFile)ccorrected->SaveAs("CorrectedPbPb.eps");
//TH1D unfoldingspectrac[fNCentralityBinAtTheEnd];
//TGraphErrors unfoldingspectracn[fNCentralityBinAtTheEnd];
TH1D *correctedspectrac = new TH1D[fNCentralityBinAtTheEnd];
TGraphErrors *correctedspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
-
-
if(fBeamType==1) {
TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
- TString titlee("corrected_centrality_bin_");
- titlee += "[";
- titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
- titlee += "_";
- titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
- titlee += "[";
- TString titleec("corrected_check_projection_bin_");
- titleec += "[";
- titleec += fLowBoundaryCentralityBinAtTheEnd[binc];
- titleec += "_";
- titleec += fHighBoundaryCentralityBinAtTheEnd[binc];
- titleec += "[";
- TString titleenameunotnormalized("Unfolded_Notnormalized_centrality_bin_");
- titleenameunotnormalized += "[";
- titleenameunotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
- titleenameunotnormalized += "_";
- titleenameunotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
- titleenameunotnormalized += "[";
- TString titleenameunormalized("Unfolded_normalized_centrality_bin_");
- titleenameunormalized += "[";
- titleenameunormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
- titleenameunormalized += "_";
- titleenameunormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
- titleenameunormalized += "[";
- TString titleenamednotnormalized("Dirrectcorrected_Notnormalized_centrality_bin_");
- titleenamednotnormalized += "[";
- titleenamednotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
- titleenamednotnormalized += "_";
- titleenamednotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
- titleenamednotnormalized += "[";
- TString titleenamednormalized("Dirrectedcorrected_normalized_centrality_bin_");
- titleenamednormalized += "[";
- titleenamednormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
- titleenamednormalized += "_";
- titleenamednormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
- titleenamednormalized += "[";
- Int_t nbEvents = 0;
- for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
- printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
- nbEvents += fNEvents[k];
- }
- Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
- Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
- printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
- Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
- Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
- printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
- cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
- cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
- TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
- ccorrectedcheck->cd(1);
- TH1D *aftersuc = (TH1D *) sparsesu->Projection(0);
- TH1D *aftersdc = (TH1D *) sparsed->Projection(0);
- aftersuc->Draw();
- aftersdc->Draw("same");
- TCanvas * ccorrectede = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
- ccorrectede->Divide(2,1);
- ccorrectede->cd(1);
- gPad->SetLogy();
- TH1D *aftersu = (TH1D *) sparsesu->Projection(1);
- CorrectFromTheWidth(aftersu);
- aftersu->SetName((const char*)titleenameunotnormalized);
- unfoldingspectrac[binc] = *aftersu;
- ccorrectede->cd(1);
- TGraphErrors* aftersun = NormalizeTH1N(aftersu,nbEvents);
- aftersun->SetTitle("");
- aftersun->GetYaxis()->SetTitleOffset(1.5);
- aftersun->GetYaxis()->SetRangeUser(0.000000001,1.0);
- aftersun->SetMarkerStyle(26);
- aftersun->SetMarkerColor(kBlue);
- aftersun->SetLineColor(kBlue);
- aftersun->Draw("AP");
- aftersun->SetName((const char*)titleenameunormalized);
- unfoldingspectracn[binc] = *aftersun;
- ccorrectede->cd(1);
- TH1D *aftersd = (TH1D *) sparsed->Projection(1);
- CorrectFromTheWidth(aftersd);
- aftersd->SetName((const char*)titleenamednotnormalized);
- correctedspectrac[binc] = *aftersd;
- ccorrectede->cd(1);
- TGraphErrors* aftersdn = NormalizeTH1N(aftersd,nbEvents);
- aftersdn->SetTitle("");
- aftersdn->GetYaxis()->SetTitleOffset(1.5);
- aftersdn->GetYaxis()->SetRangeUser(0.000000001,1.0);
- aftersdn->SetMarkerStyle(25);
- aftersdn->SetMarkerColor(kBlack);
- aftersdn->SetLineColor(kBlack);
- aftersdn->Draw("P");
- aftersdn->SetName((const char*)titleenamednormalized);
- correctedspectracn[binc] = *aftersdn;
- TLegend *legcorrectedud = new TLegend(0.4,0.6,0.89,0.89);
- legcorrectedud->AddEntry(aftersun,"Corrected","p");
- legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
- legcorrectedud->Draw("same");
- ccorrectedallspectra->cd(1);
- gPad->SetLogy();
- TH1D *aftersunn = (TH1D *) aftersun->Clone();
- aftersunn->SetMarkerStyle(stylee[binc]);
- aftersunn->SetMarkerColor(colorr[binc]);
- if(binc==0) aftersunn->Draw("AP");
- else aftersunn->Draw("P");
- legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
- ccorrectedallspectra->cd(2);
- gPad->SetLogy();
- TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
- aftersdnn->SetMarkerStyle(stylee[binc]);
- aftersdnn->SetMarkerColor(colorr[binc]);
- if(binc==0) aftersdnn->Draw("AP");
- else aftersdnn->Draw("P");
- legtotalg->AddEntry(aftersdnn,(const char*) titlee,"p");
- ccorrectede->cd(2);
- TH1D* ratiocorrectedbinc = (TH1D*)aftersu->Clone();
- TString titleee("ratiocorrected_bin_");
- titleee += binc;
- ratiocorrectedbinc->SetName((const char*) titleee);
- ratiocorrectedbinc->SetTitle("");
- ratiocorrectedbinc->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
- ratiocorrectedbinc->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- ratiocorrectedbinc->Divide(aftersu,aftersd,1,1);
- ratiocorrectedbinc->SetStats(0);
- ratiocorrectedbinc->Draw();
+ TString titlee("corrected_centrality_bin_");
+ titlee += "[";
+ titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
+ titlee += "_";
+ titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
+ titlee += "[";
+ TString titleec("corrected_check_projection_bin_");
+ titleec += "[";
+ titleec += fLowBoundaryCentralityBinAtTheEnd[binc];
+ titleec += "_";
+ titleec += fHighBoundaryCentralityBinAtTheEnd[binc];
+ titleec += "[";
+ TString titleenameunotnormalized("Unfolded_Notnormalized_centrality_bin_");
+ titleenameunotnormalized += "[";
+ titleenameunotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
+ titleenameunotnormalized += "_";
+ titleenameunotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
+ titleenameunotnormalized += "[";
+ TString titleenameunormalized("Unfolded_normalized_centrality_bin_");
+ titleenameunormalized += "[";
+ titleenameunormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
+ titleenameunormalized += "_";
+ titleenameunormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
+ titleenameunormalized += "[";
+ TString titleenamednotnormalized("Dirrectcorrected_Notnormalized_centrality_bin_");
+ titleenamednotnormalized += "[";
+ titleenamednotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
+ titleenamednotnormalized += "_";
+ titleenamednotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
+ titleenamednotnormalized += "[";
+ TString titleenamednormalized("Dirrectedcorrected_normalized_centrality_bin_");
+ titleenamednormalized += "[";
+ titleenamednormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
+ titleenamednormalized += "_";
+ titleenamednormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
+ titleenamednormalized += "[";
+ Int_t nbEvents = 0;
+ for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
+ printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
+ nbEvents += fNEvents[k];
+ }
+ Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
+ Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
+ printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
+ Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
+ Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
+ printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
+ cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
+ cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
+ TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
+ ccorrectedcheck->cd(1);
+ TH1D *aftersuc = (TH1D *) sparsesu->Projection(0);
+ TH1D *aftersdc = (TH1D *) sparsed->Projection(0);
+ aftersuc->Draw();
+ aftersdc->Draw("same");
+ TCanvas * ccorrectede = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
+ ccorrectede->Divide(2,1);
+ ccorrectede->cd(1);
+ gPad->SetLogy();
+ TH1D *aftersu = (TH1D *) sparsesu->Projection(1);
+ CorrectFromTheWidth(aftersu);
+ aftersu->SetName((const char*)titleenameunotnormalized);
+ unfoldingspectrac[binc] = *aftersu;
+ ccorrectede->cd(1);
+ TGraphErrors* aftersun = NormalizeTH1N(aftersu,nbEvents);
+ aftersun->SetTitle("");
+ aftersun->GetYaxis()->SetTitleOffset(1.5);
+ aftersun->GetYaxis()->SetRangeUser(0.000000001,1.0);
+ aftersun->SetMarkerStyle(26);
+ aftersun->SetMarkerColor(kBlue);
+ aftersun->SetLineColor(kBlue);
+ aftersun->Draw("AP");
+ aftersun->SetName((const char*)titleenameunormalized);
+ unfoldingspectracn[binc] = *aftersun;
+ ccorrectede->cd(1);
+ TH1D *aftersd = (TH1D *) sparsed->Projection(1);
+ CorrectFromTheWidth(aftersd);
+ aftersd->SetName((const char*)titleenamednotnormalized);
+ correctedspectrac[binc] = *aftersd;
+ ccorrectede->cd(1);
+ TGraphErrors* aftersdn = NormalizeTH1N(aftersd,nbEvents);
+ aftersdn->SetTitle("");
+ aftersdn->GetYaxis()->SetTitleOffset(1.5);
+ aftersdn->GetYaxis()->SetRangeUser(0.000000001,1.0);
+ aftersdn->SetMarkerStyle(25);
+ aftersdn->SetMarkerColor(kBlack);
+ aftersdn->SetLineColor(kBlack);
+ aftersdn->Draw("P");
+ aftersdn->SetName((const char*)titleenamednormalized);
+ correctedspectracn[binc] = *aftersdn;
+ TLegend *legcorrectedud = new TLegend(0.4,0.6,0.89,0.89);
+ legcorrectedud->AddEntry(aftersun,"Corrected","p");
+ legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
+ legcorrectedud->Draw("same");
+ ccorrectedallspectra->cd(1);
+ gPad->SetLogy();
+ TH1D *aftersunn = (TH1D *) aftersun->Clone();
+ aftersunn->SetMarkerStyle(stylee[binc]);
+ aftersunn->SetMarkerColor(colorr[binc]);
+ if(binc==0) aftersunn->Draw("AP");
+ else aftersunn->Draw("P");
+ legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
+ ccorrectedallspectra->cd(2);
+ gPad->SetLogy();
+ TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
+ aftersdnn->SetMarkerStyle(stylee[binc]);
+ aftersdnn->SetMarkerColor(colorr[binc]);
+ if(binc==0) aftersdnn->Draw("AP");
+ else aftersdnn->Draw("P");
+ legtotalg->AddEntry(aftersdnn,(const char*) titlee,"p");
+ ccorrectede->cd(2);
+ TH1D* ratiocorrectedbinc = (TH1D*)aftersu->Clone();
+ TString titleee("ratiocorrected_bin_");
+ titleee += binc;
+ ratiocorrectedbinc->SetName((const char*) titleee);
+ ratiocorrectedbinc->SetTitle("");
+ ratiocorrectedbinc->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
+ ratiocorrectedbinc->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+ ratiocorrectedbinc->Divide(aftersu,aftersd,1,1);
+ ratiocorrectedbinc->SetStats(0);
+ ratiocorrectedbinc->Draw();
}
ccorrectedallspectra->cd(1);
cenaxisa->SetRange(0,nbbin);
cenaxisb->SetRange(0,nbbin);
-
+ if(fWriteToFile) ccorrectedallspectra->SaveAs("CorrectedPbPb.eps");
}
// Dump to file if needed
alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
alltogetherCorrection->Write();
for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
- unfoldingspectrac[binc].Write();
- unfoldingspectracn[binc].Write();
- correctedspectrac[binc].Write();
- correctedspectracn[binc].Write();
+ unfoldingspectrac[binc].Write();
+ unfoldingspectracn[binc].Write();
+ correctedspectrac[binc].Write();
+ correctedspectracn[binc].Write();
}
out->Close(); delete out;
}
}
-
-
-
return kTRUE;
}
// Subtract hadron background
/////////////////////////////////
AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
+ AliCFDataGrid *unnormalizedRawSpectrum = 0x0;
+ TGraphErrors *gNormalizedRawSpectrum = 0x0;
if(subtractcontamination) {
- dataspectrumaftersubstraction = SubtractBackground(kTRUE);
- dataGridAfterFirstSteps = dataspectrumaftersubstraction;
+ if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE);
+ else dataspectrumaftersubstraction = GetRawBspectra2ndMethod();
+ unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
+ dataGridAfterFirstSteps = dataspectrumaftersubstraction;
+ gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
}
+ printf("after normalize getting IP \n");
+
/////////////////////////////////////////////////////////////////////////////////////////
// Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
/////////////////////////////////////////////////////////////////////////////////////////
if(fEfficiencyFunction){
dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
- dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
+ dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
}
else if(dataContainerV0){
dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
}
+
///////////////
// Unfold
//////////////
/////////////////////
// Simply correct
////////////////////
+
AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
//////////
// Plot
//////////
+
if(fDebugLevel > 0.0) {
+
+ Int_t ptpr = 0;
+ if(fBeamType==0) ptpr=0;
+ if(fBeamType==1) ptpr=1;
TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
ccorrected->Divide(2,1);
legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
legcorrected->Draw("same");
ccorrected->cd(2);
- TH1D *correctedTH1D = correctedspectrum->Projection(0);
- TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(0);
+ TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
+ TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
ratiocorrected->SetName("ratiocorrected");
ratiocorrected->SetTitle("");
ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
ratiocorrected->SetStats(0);
ratiocorrected->Draw();
+ if(fWriteToFile) ccorrected->SaveAs("CorrectedBeauty.eps");
+ if(fBeamType == 0){
+ if(fNonHFEsyst){
+ CalculateNonHFEsyst(0);
+ }
+ }
// Dump to file if needed
if(fDumpToFile) {
-
+ // to do centrality dependent
+
TFile *out;
- if(fNonHFEmode == 1){
- out = new TFile("finalSpectrumLow.root","recreate");
- }
- else if(fNonHFEmode == 2){
- out = new TFile("finalSpectrumUp.root","recreate");
- }
- else{
- out = new TFile("finalSpectrum.root","recreate");
- }
+ out = new TFile("finalSpectrum.root","recreate");
out->cd();
//
correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
alltogetherCorrection->Write();
//
+ if(unnormalizedRawSpectrum) {
+ unnormalizedRawSpectrum->SetName("beautyAfterIP");
+ unnormalizedRawSpectrum->Write();
+ }
+
+ if(gNormalizedRawSpectrum){
+ gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
+ gNormalizedRawSpectrum->Write();
+ }
+
+ if(fBeamType==0) {
+ Int_t countpp=0;
+ fEfficiencyCharmSigD[countpp]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",countpp));
+ fEfficiencyCharmSigD[countpp]->SetName(Form("IPEfficiencyForCharmSigCent%i",countpp));
+ fEfficiencyCharmSigD[countpp]->Write();
+ fEfficiencyBeautySigD[countpp]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",countpp));
+ fEfficiencyBeautySigD[countpp]->SetName(Form("IPEfficiencyForBeautySigCent%i",countpp));
+ fEfficiencyBeautySigD[countpp]->Write();
+ fCharmEff[countpp]->SetTitle(Form("IPEfficiencyForCharmCent%i",countpp));
+ fCharmEff[countpp]->SetName(Form("IPEfficiencyForCharmCent%i",countpp));
+ fCharmEff[countpp]->Write();
+ fBeautyEff[countpp]->SetTitle(Form("IPEfficiencyForBeautyCent%i",countpp));
+ fBeautyEff[countpp]->SetName(Form("IPEfficiencyForBeautyCent%i",countpp));
+ fBeautyEff[countpp]->Write();
+ fConversionEff[countpp]->SetTitle(Form("IPEfficiencyForConversionCent%i",countpp));
+ fConversionEff[countpp]->SetName(Form("IPEfficiencyForConversionCent%i",countpp));
+ fConversionEff[countpp]->Write();
+ fNonHFEEff[countpp]->SetTitle(Form("IPEfficiencyForNonHFECent%i",countpp));
+ fNonHFEEff[countpp]->SetName(Form("IPEfficiencyForNonHFECent%i",countpp));
+ fNonHFEEff[countpp]->Write();
+ }
+
+ if(fBeamType==1) {
+
+ TGraphErrors* correctedspectrumDc[kCentrality];
+ TGraphErrors* alltogetherspectrumDc[kCentrality];
+ for(Int_t i=0;i<kCentrality-2;i++)
+ {
+ correctedspectrum->GetAxis(0)->SetRange(i+1,i+1);
+ correctedspectrumDc[i] = Normalize(correctedspectrum,i);
+ if(correctedspectrumDc[i]){
+ correctedspectrumDc[i]->SetTitle(Form("UnfoldingCorrectedSpectrum_%i",i));
+ correctedspectrumDc[i]->SetName(Form("UnfoldingCorrectedSpectrum_%i",i));
+ correctedspectrumDc[i]->Write();
+ }
+ alltogetherCorrection->GetAxis(0)->SetRange(i+1,i+1);
+ alltogetherspectrumDc[i] = Normalize(alltogetherCorrection,i);
+ if(alltogetherspectrumDc[i]){
+ alltogetherspectrumDc[i]->SetTitle(Form("AlltogetherSpectrum_%i",i));
+ alltogetherspectrumDc[i]->SetName(Form("AlltogetherSpectrum_%i",i));
+ alltogetherspectrumDc[i]->Write();
+ }
+
+ TH1D *centrcrosscheck = correctedspectrum->Projection(0);
+ centrcrosscheck->SetTitle(Form("centrality_%i",i));
+ centrcrosscheck->SetName(Form("centrality_%i",i));
+ centrcrosscheck->Write();
+
+ TH1D *correctedTH1Dc = correctedspectrum->Projection(ptpr);
+ TH1D *alltogetherTH1Dc = (TH1D *) alltogetherCorrection->Project(ptpr);
+
+ TH1D *centrcrosscheck2 = (TH1D *) alltogetherCorrection->Project(0);
+ centrcrosscheck2->SetTitle(Form("centrality2_%i",i));
+ centrcrosscheck2->SetName(Form("centrality2_%i",i));
+ centrcrosscheck2->Write();
+
+ TH1D* ratiocorrectedc = (TH1D*)correctedTH1D->Clone();
+ ratiocorrectedc->Divide(correctedTH1Dc,alltogetherTH1Dc,1,1);
+ ratiocorrectedc->SetTitle(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
+ ratiocorrectedc->SetName(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
+ ratiocorrectedc->Write();
+
+ fEfficiencyCharmSigD[i]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",i));
+ fEfficiencyCharmSigD[i]->SetName(Form("IPEfficiencyForCharmSigCent%i",i));
+ fEfficiencyCharmSigD[i]->Write();
+ fEfficiencyBeautySigD[i]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",i));
+ fEfficiencyBeautySigD[i]->SetName(Form("IPEfficiencyForBeautySigCent%i",i));
+ fEfficiencyBeautySigD[i]->Write();
+ fCharmEff[i]->SetTitle(Form("IPEfficiencyForCharmCent%i",i));
+ fCharmEff[i]->SetName(Form("IPEfficiencyForCharmCent%i",i));
+ fCharmEff[i]->Write();
+ fBeautyEff[i]->SetTitle(Form("IPEfficiencyForBeautyCent%i",i));
+ fBeautyEff[i]->SetName(Form("IPEfficiencyForBeautyCent%i",i));
+ fBeautyEff[i]->Write();
+ fConversionEff[i]->SetTitle(Form("IPEfficiencyForConversionCent%i",i));
+ fConversionEff[i]->SetName(Form("IPEfficiencyForConversionCent%i",i));
+ fConversionEff[i]->Write();
+ fNonHFEEff[i]->SetTitle(Form("IPEfficiencyForNonHFECent%i",i));
+ fNonHFEEff[i]->SetName(Form("IPEfficiencyForNonHFECent%i",i));
+ fNonHFEEff[i]->Write();
+ }
+
+ }
+
out->Close();
delete out;
}
-
-
}
return kTRUE;
//
// Apply background subtraction
//
-
+
+ Int_t ptpr = 0;
+ Int_t nbins=1;
+ if(fBeamType==0)
+ {
+ ptpr=0;
+ nbins=1;
+ }
+ if(fBeamType==1)
+ {
+ ptpr=1;
+ nbins=2;
+ }
+
// Raw spectrum
AliCFContainer *dataContainer = GetContainer(kDataContainer);
if(!dataContainer){
AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
+
// Background Estimate
AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
if(!backgroundContainer){
if(!fInclusiveSpectrum){
//Background subtraction for IP analysis
- TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(0);
+
+ TH1D *incElecCent[kCentrality-1];
+ TH1D *charmCent[kCentrality-1];
+ TH1D *convCent[kCentrality-1];
+ TH1D *nonHFECent[kCentrality-1];
+ TH1D *subtractedCent[kCentrality-1];
+ TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
CorrectFromTheWidth(measuredTH1Draw);
- TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",600,500);
+ if(fBeamType==1){
+ THnSparseF* sparseIncElec = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
+ for(Int_t icent = 1; icent < kCentrality-1; icent++){
+ sparseIncElec->GetAxis(0)->SetRange(icent,icent);
+ incElecCent[icent-1] = (TH1D *) sparseIncElec->Projection(ptpr);
+ CorrectFromTheWidth(incElecCent[icent-1]);
+ }
+ }
+ TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",500,400);
rawspectra->cd();
+ rawspectra->SetLogy();
+ gStyle->SetOptStat(0);
+ TLegend *lRaw = new TLegend(0.55,0.55,0.85,0.85);
measuredTH1Draw->SetMarkerStyle(20);
measuredTH1Draw->Draw();
+ measuredTH1Draw->GetXaxis()->SetRangeUser(0.0,7.9);
+ lRaw->AddEntry(measuredTH1Draw,"measured raw spectrum");
+ TH1D* htemp;
+ Int_t* bins=new Int_t[2];
if(fIPanaHadronBgSubtract){
// Hadron background
- printf("Hadron background for IP analysis subtracted!\n");
- Int_t* bins=new Int_t[1];
- TH1D* htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
- bins[0]=htemp->GetNbinsX();
- AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",1,bins);
+ printf("Hadron background for IP analysis subtracted!\n");
+ if(fBeamType==0)
+ {
+
+ htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
+ bins[0]=htemp->GetNbinsX();
+ }
+ if(fBeamType==1)
+ {
+ htemp = (TH1D *) fHadronEffbyIPcut->Projection(0);
+ bins[0]=htemp->GetNbinsX();
+ htemp = (TH1D *) fHadronEffbyIPcut->Projection(1);
+ bins[1]=htemp->GetNbinsX();
+ }
+ AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins);
hbgContainer->SetGrid(fHadronEffbyIPcut);
backgroundGrid->Multiply(hbgContainer,1);
// draw raw hadron bg spectra
- TH1D *hadronbg= (TH1D *) backgroundGrid->Project(0);
+ TH1D *hadronbg= (TH1D *) backgroundGrid->Project(ptpr);
CorrectFromTheWidth(hadronbg);
hadronbg->SetMarkerColor(7);
hadronbg->SetMarkerStyle(20);
rawspectra->cd();
hadronbg->Draw("samep");
+ lRaw->AddEntry(hadronbg,"hadrons");
// subtract hadron contamination
spectrumSubtracted->Add(backgroundGrid,-1.0);
}
printf("Charm background for IP analysis subtracted!\n");
AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
// draw charm bg spectra
- TH1D *charmbg= (TH1D *) charmbgContainer->Project(0);
+ TH1D *charmbg= (TH1D *) charmbgContainer->Project(ptpr);
CorrectFromTheWidth(charmbg);
charmbg->SetMarkerColor(3);
charmbg->SetMarkerStyle(20);
rawspectra->cd();
charmbg->Draw("samep");
+ lRaw->AddEntry(charmbg,"charm elecs");
// subtract charm background
spectrumSubtracted->Add(charmbgContainer,-1.0);
+ if(fBeamType==1){
+ THnSparseF* sparseCharmElec = (THnSparseF *) charmbgContainer->GetGrid();
+ for(Int_t icent = 1; icent < kCentrality-1; icent++){
+ sparseCharmElec->GetAxis(0)->SetRange(icent,icent);
+ charmCent[icent-1] = (TH1D *) sparseCharmElec->Projection(ptpr);
+ CorrectFromTheWidth(charmCent[icent-1]);
+ }
+ }
}
if(fIPanaConversionBgSubtract){
// Conversion background
- AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
+ AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
// draw conversion bg spectra
- TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(0);
+ TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(ptpr);
CorrectFromTheWidth(conversionbg);
conversionbg->SetMarkerColor(4);
conversionbg->SetMarkerStyle(20);
rawspectra->cd();
conversionbg->Draw("samep");
+ lRaw->AddEntry(conversionbg,"conversion elecs");
// subtract conversion background
spectrumSubtracted->Add(conversionbgContainer,-1.0);
- printf("Conversion background subtraction is preliminary!\n");
+ if(fBeamType==1){
+ THnSparseF* sparseconvElec = (THnSparseF *) conversionbgContainer->GetGrid();
+ for(Int_t icent = 1; icent < kCentrality-1; icent++){
+ sparseconvElec->GetAxis(0)->SetRange(icent,icent);
+ convCent[icent-1] = (TH1D *) sparseconvElec->Projection(ptpr);
+ CorrectFromTheWidth(convCent[icent-1]);
+ }
+ }
}
if(fIPanaNonHFEBgSubtract){
// NonHFE background
AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
// draw Dalitz/dielectron bg spectra
- TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(0);
+ TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(ptpr);
CorrectFromTheWidth(nonhfebg);
nonhfebg->SetMarkerColor(6);
nonhfebg->SetMarkerStyle(20);
rawspectra->cd();
nonhfebg->Draw("samep");
+ lRaw->AddEntry(nonhfebg,"non-HF elecs");
// subtract Dalitz/dielectron background
spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
- printf("Non HFE background subtraction is preliminary!\n");
+ if(fBeamType==1){
+ THnSparseF* sparseNonHFEElec = (THnSparseF *) nonHFEbgContainer->GetGrid();
+ for(Int_t icent = 1; icent < kCentrality-1; icent++){
+ sparseNonHFEElec->GetAxis(0)->SetRange(icent,icent);
+ nonHFECent[icent-1] = (TH1D *) sparseNonHFEElec->Projection(ptpr);
+ CorrectFromTheWidth(nonHFECent[icent-1]);
+ }
+ }
}
- TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(0);
+
+ TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(ptpr);
CorrectFromTheWidth(rawbgsubtracted);
rawbgsubtracted->SetMarkerStyle(24);
rawspectra->cd();
+ lRaw->AddEntry(rawbgsubtracted,"subtracted raw spectrum");
rawbgsubtracted->Draw("samep");
+ lRaw->Draw("SAME");
+ gPad->SetGrid();
+ //rawspectra->SaveAs("rawspectra.eps");
+
+ if(fBeamType==1){
+ THnSparseF* sparseSubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
+ for(Int_t icent = 1; icent < kCentrality-1; icent++){
+ sparseSubtracted->GetAxis(0)->SetRange(icent,icent);
+ subtractedCent[icent-1] = (TH1D *) sparseSubtracted->Projection(ptpr);
+ CorrectFromTheWidth(subtractedCent[icent-1]);
+ }
+
+ TLegend *lCentRaw = new TLegend(0.55,0.55,0.85,0.85);
+ TCanvas *centRaw = new TCanvas("centRaw","centRaw",1000,800);
+ centRaw->Divide(3,3);
+ for(Int_t icent = 1; icent < kCentrality-1; icent++){
+ centRaw->cd(icent);
+ gPad->SetLogx();
+ gPad->SetLogy();
+ incElecCent[icent-1]->GetXaxis()->SetRangeUser(0.4,8.);
+ incElecCent[icent-1]->Draw("p");
+ incElecCent[icent-1]->SetMarkerColor(1);
+ incElecCent[icent-1]->SetMarkerStyle(20);
+ charmCent[icent-1]->Draw("samep");
+ charmCent[icent-1]->SetMarkerColor(3);
+ charmCent[icent-1]->SetMarkerStyle(20);
+ convCent[icent-1]->Draw("samep");
+ convCent[icent-1]->SetMarkerColor(4);
+ convCent[icent-1]->SetMarkerStyle(20);
+ nonHFECent[icent-1]->Draw("samep");
+ nonHFECent[icent-1]->SetMarkerColor(6);
+ nonHFECent[icent-1]->SetMarkerStyle(20);
+ subtractedCent[icent-1]->Draw("samep");
+ subtractedCent[icent-1]->SetMarkerStyle(24);
+ if(icent == 1){
+ lCentRaw->AddEntry(incElecCent[0],"inclusive electron spectrum");
+ lCentRaw->AddEntry(charmCent[0],"charm elecs");
+ lCentRaw->AddEntry(convCent[0],"conversion elecs");
+ lCentRaw->AddEntry(nonHFECent[0],"non-HF elecs");
+ lCentRaw->AddEntry(subtractedCent[0],"subtracted electron spectrum");
+ lCentRaw->Draw("SAME");
+ }
+ }
+ }
+
+ delete[] bins;
+
}
else{
// Subtract
if(fDebugLevel > 0) {
- Int_t ptpr;
- if(fBeamType==0) ptpr=0;
- if(fBeamType==1) ptpr=1;
+ Int_t ptprd;
+ if(fBeamType==0) ptprd=0;
+ if(fBeamType==1) ptprd=1;
TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
cbackgroundsubtraction->Divide(3,1);
cbackgroundsubtraction->cd(1);
gPad->SetLogy();
- TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptpr);
- TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
+ TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptprd);
+ TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptprd);
CorrectFromTheWidth(measuredTH1Daftersubstraction);
CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
measuredTH1Daftersubstraction->SetStats(0);
}
ratiomeasuredcontamination->Draw("P");
cbackgroundsubtraction->cd(3);
- TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptpr);
+ TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptprd);
CorrectFromTheWidth(measuredTH1background);
measuredTH1background->SetStats(0);
measuredTH1background->SetTitle("");
measuredTH1background->SetMarkerColor(kBlack);
measuredTH1background->SetLineColor(kBlack);
measuredTH1background->Draw();
+ if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.eps");
if(fBeamType==1) {
Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
for(Int_t binc = 0; binc < nbbin; binc++){
- TString titlee("BackgroundSubtraction_centrality_bin_");
- titlee += binc;
- TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
- cbackground->Divide(2,1);
- cbackground->cd(1);
- gPad->SetLogy();
- cenaxisa->SetRange(binc+1,binc+1);
- cenaxisb->SetRange(binc+1,binc+1);
- TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
- TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
- CorrectFromTheWidth(aftersubstraction);
- CorrectFromTheWidth(beforesubstraction);
- aftersubstraction->SetStats(0);
- aftersubstraction->SetTitle((const char*)titlee);
- aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
- aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- aftersubstraction->SetMarkerStyle(25);
- aftersubstraction->SetMarkerColor(kBlack);
- aftersubstraction->SetLineColor(kBlack);
- beforesubstraction->SetStats(0);
- beforesubstraction->SetTitle((const char*)titlee);
- beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
- beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- beforesubstraction->SetMarkerStyle(24);
- beforesubstraction->SetMarkerColor(kBlue);
- beforesubstraction->SetLineColor(kBlue);
- aftersubstraction->Draw();
- beforesubstraction->Draw("same");
- TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
- lega->AddEntry(beforesubstraction,"With hadron contamination","p");
- lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
- lega->Draw("same");
- cbackgrounde->cd(1);
- gPad->SetLogy();
- TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
- aftersubtractionn->SetMarkerStyle(stylee[binc]);
- aftersubtractionn->SetMarkerColor(colorr[binc]);
- if(binc==0) aftersubtractionn->Draw();
- else aftersubtractionn->Draw("same");
- legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
- cbackgrounde->cd(2);
- gPad->SetLogy();
- TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
- aftersubtractionng->SetMarkerStyle(stylee[binc]);
- aftersubtractionng->SetMarkerColor(colorr[binc]);
- if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
- if(binc==0) aftersubtractionng->Draw();
- else aftersubtractionng->Draw("same");
- legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
- cbackground->cd(2);
- TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
- ratiocontamination->SetName("ratiocontamination");
- ratiocontamination->SetTitle((const char*)titlee);
- ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
- ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- ratiocontamination->Add(aftersubstraction,-1.0);
- ratiocontamination->Divide(beforesubstraction);
- Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
- for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
- ratiocontamination->SetBinError(nbinpt+1,0.0);
- }
- ratiocontamination->SetStats(0);
- ratiocontamination->SetMarkerStyle(26);
- ratiocontamination->SetMarkerColor(kBlack);
- ratiocontamination->SetLineColor(kBlack);
- ratiocontamination->Draw("P");
-
+ TString titlee("BackgroundSubtraction_centrality_bin_");
+ titlee += binc;
+ TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
+ cbackground->Divide(2,1);
+ cbackground->cd(1);
+ gPad->SetLogy();
+ cenaxisa->SetRange(binc+1,binc+1);
+ cenaxisb->SetRange(binc+1,binc+1);
+ TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
+ TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
+ CorrectFromTheWidth(aftersubstraction);
+ CorrectFromTheWidth(beforesubstraction);
+ aftersubstraction->SetStats(0);
+ aftersubstraction->SetTitle((const char*)titlee);
+ aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+ aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+ aftersubstraction->SetMarkerStyle(25);
+ aftersubstraction->SetMarkerColor(kBlack);
+ aftersubstraction->SetLineColor(kBlack);
+ beforesubstraction->SetStats(0);
+ beforesubstraction->SetTitle((const char*)titlee);
+ beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+ beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+ beforesubstraction->SetMarkerStyle(24);
+ beforesubstraction->SetMarkerColor(kBlue);
+ beforesubstraction->SetLineColor(kBlue);
+ aftersubstraction->Draw();
+ beforesubstraction->Draw("same");
+ TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
+ lega->AddEntry(beforesubstraction,"With hadron contamination","p");
+ lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
+ lega->Draw("same");
+ cbackgrounde->cd(1);
+ gPad->SetLogy();
+ TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
+ aftersubtractionn->SetMarkerStyle(stylee[binc]);
+ aftersubtractionn->SetMarkerColor(colorr[binc]);
+ if(binc==0) aftersubtractionn->Draw();
+ else aftersubtractionn->Draw("same");
+ legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
+ cbackgrounde->cd(2);
+ gPad->SetLogy();
+ TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
+ aftersubtractionng->SetMarkerStyle(stylee[binc]);
+ aftersubtractionng->SetMarkerColor(colorr[binc]);
+ if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
+ if(binc==0) aftersubtractionng->Draw();
+ else aftersubtractionng->Draw("same");
+ legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
+ cbackground->cd(2);
+ TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
+ ratiocontamination->SetName("ratiocontamination");
+ ratiocontamination->SetTitle((const char*)titlee);
+ ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
+ ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+ ratiocontamination->Add(aftersubstraction,-1.0);
+ ratiocontamination->Divide(beforesubstraction);
+ Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
+ for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
+ ratiocontamination->SetBinError(nbinpt+1,0.0);
+ }
+ ratiocontamination->SetStats(0);
+ ratiocontamination->SetMarkerStyle(26);
+ ratiocontamination->SetMarkerColor(kBlack);
+ ratiocontamination->SetLineColor(kBlack);
+ ratiocontamination->Draw("P");
}
cbackgrounde->cd(1);
cenaxisa->SetRange(0,nbbin);
cenaxisb->SetRange(0,nbbin);
-
+ if(fWriteToFile) cbackgrounde->SaveAs("BackgroundSubtractedPbPb.eps");
}
-
-
}
return spectrumSubtracted;
//
// calculate charm background
//
+ Int_t ptpr = 0;
+ Int_t nDim = 1;
+ if(fBeamType==0)
+ {
+ ptpr=0;
+ }
+ if(fBeamType==1)
+ {
+ ptpr=1;
+ nDim=2;
+ }
Double_t evtnorm=0;
- if(fNMCbgEvents) evtnorm= double(fNEvents[0])/double(fNMCbgEvents);
+ if(fNMCbgEvents[0]) evtnorm= double(fNEvents[0])/double(fNMCbgEvents[0]);
AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
if(!mcContainer){
AliCFDataGrid *charmBackgroundGrid= 0x0;
charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
- TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
- Int_t* bins=new Int_t[1];
+ TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
+ Int_t* bins=new Int_t[2];
bins[0]=charmbgaftertofpid->GetNbinsX();
- AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",1,bins);
+
+ if(fBeamType==1)
+ {
+ bins[0]=12;
+ charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(1);
+ bins[1]=charmbgaftertofpid->GetNbinsX();
+
+ }
+
+ AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
- TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0);
+ TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(ptpr);
+
+ charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.);
+
+ Int_t *nBinpp=new Int_t[1];
+ Int_t* binspp=new Int_t[1];
+ binspp[0]=charmbgaftertofpid->GetNbinsX(); // number of pt bins
+
+ Int_t *nBinPbPb=new Int_t[2];
+ Int_t* binsPbPb=new Int_t[2];
+ binsPbPb[1]=charmbgaftertofpid->GetNbinsX(); // number of pt bins
+ binsPbPb[0]=12;
- charmBackgroundGrid->Multiply(ipWeightedCharmContainer,evtnorm);
- TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0);
+ Int_t looppt=binspp[0];
+ if(fBeamType==1) looppt=binsPbPb[1];
- AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",1,bins);
+ for(Long_t iBin=1; iBin<= looppt;iBin++){
+ if(fBeamType==0)
+ {
+ nBinpp[0]=iBin;
+ charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
+ charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
+ }
+ if(fBeamType==1)
+ {
+ // loop over centrality
+ for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
+ nBinPbPb[0]=iiBin;
+ nBinPbPb[1]=iBin;
+ Double_t evtnormPbPb=0;
+ if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
+ charmBackgroundGrid->SetElementError(nBinPbPb, charmBackgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
+ charmBackgroundGrid->SetElement(nBinPbPb,charmBackgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
+ }
+ }
+ }
+
+ TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(ptpr);
+
+ AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors
- TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0);
+ TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(ptpr);
charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
- TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0);
+ TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(ptpr);
// Efficiency (set efficiency to 1 for only folding)
AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
efficiencyD->CalculateEfficiency(0,0);
- // Folding
- Int_t nDim = 1;
+ // Folding
+ if(fBeamType==0)nDim = 1;
+ if(fBeamType==1)nDim = 2;
AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
folding.SetMaxNumberOfIterations(1);
folding.Unfold();
THnSparse* result1= folding.GetEstMeasured(); // folded spectra
THnSparse* result=(THnSparse*)result1->Clone();
charmBackgroundGrid->SetGrid(result);
- TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0);
+ TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(ptpr);
//Charm background evaluation plots
cCharmBgEval->Divide(3,1);
cCharmBgEval->cd(1);
- charmbgaftertofpid->Scale(evtnorm);
+
+ if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
+ if(fBeamType==1)
+ {
+ Double_t evtnormPbPb=0;
+ for(Int_t kCentr=0;kCentr<bins[0];kCentr++)
+ {
+ if(fNMCbgEvents[kCentr]) evtnormPbPb= evtnormPbPb+double(fNEvents[kCentr])/double(fNMCbgEvents[kCentr]);
+ }
+ charmbgaftertofpid->Scale(evtnormPbPb);
+ }
+
CorrectFromTheWidth(charmbgaftertofpid);
charmbgaftertofpid->SetMarkerStyle(25);
charmbgaftertofpid->Draw("p");
legcharmbg2->Draw("same");
CorrectStatErr(charmBackgroundGrid);
+ if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
+
+ delete[] bins;
+ delete[] nBinpp;
+ delete[] binspp;
+ delete[] nBinPbPb;
+ delete[] binsPbPb;
+
+ if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
return charmBackgroundGrid;
}
//
// calculate conversion background
//
-
+
Double_t evtnorm[1] = {0.0};
- if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
+ if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
printf("check event!!! %lf \n",evtnorm[0]);
-
- // Background Estimate
- AliCFContainer *backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
+
+ AliCFContainer *backgroundContainer = 0x0;
+
+ if(fNonHFEsyst){
+ backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0][0]->Clone();
+ for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
+ backgroundContainer->Add(fConvSourceContainer[iSource][0][0]); // make centrality dependent
+ }
+ }
+ else{
+ // Background Estimate
+ backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
+ }
if(!backgroundContainer){
AliError("MC background container not found");
return NULL;
}
-
+
Int_t stepbackground = 3;
+
AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
- backgroundGrid->Scale(evtnorm);
+ Int_t *nBinpp=new Int_t[1];
+ Int_t* binspp=new Int_t[1];
+ binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins
+
+ Int_t *nBinPbPb=new Int_t[2];
+ Int_t* binsPbPb=new Int_t[2];
+ binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins
+ binsPbPb[0]=12;
+
+ Int_t looppt=binspp[0];
+ if(fBeamType==1) looppt=binsPbPb[1];
+
+ for(Long_t iBin=1; iBin<= looppt;iBin++){
+ if(fBeamType==0)
+ {
+ nBinpp[0]=iBin;
+ backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
+ backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
+ }
+ if(fBeamType==1)
+ {
+ // loop over centrality
+ for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
+ nBinPbPb[0]=iiBin;
+ nBinPbPb[1]=iBin;
+ Double_t evtnormPbPb=0;
+ if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
+ backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
+ backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
+ }
+ }
+ }
+ //end of workaround for statistical errors
- Int_t* bins=new Int_t[1];
- bins[0]=fConversionEff->GetNbinsX();
- AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
+ AliCFDataGrid *weightedConversionContainer;
+ if(fBeamType==0) weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,binspp);
+ else weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",2,binsPbPb);
weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
backgroundGrid->Multiply(weightedConversionContainer,1.0);
-
- CorrectStatErr(backgroundGrid);
-
+
+ delete[] nBinpp;
+ delete[] binspp;
+ delete[] nBinPbPb;
+ delete[] binsPbPb;
+
return backgroundGrid;
}
//
Double_t evtnorm[1] = {0.0};
- if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
+ if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
printf("check event!!! %lf \n",evtnorm[0]);
-
- // Background Estimate
- AliCFContainer *backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
+
+ AliCFContainer *backgroundContainer = 0x0;
+ if(fNonHFEsyst){
+ backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0][0]->Clone();
+ for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
+ backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0]);
+ }
+ }
+ else{
+ // Background Estimate
+ backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
+ }
if(!backgroundContainer){
AliError("MC background container not found");
return NULL;
}
-
-
+
+
Int_t stepbackground = 3;
+
AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
- backgroundGrid->Scale(evtnorm);
+ Int_t *nBinpp=new Int_t[1];
+ Int_t* binspp=new Int_t[1];
+ binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins
- Int_t* bins=new Int_t[1];
- bins[0]=fNonHFEEff->GetNbinsX();
- AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
+ Int_t *nBinPbPb=new Int_t[2];
+ Int_t* binsPbPb=new Int_t[2];
+ binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins
+ binsPbPb[0]=12;
+
+ Int_t looppt=binspp[0];
+ if(fBeamType==1) looppt=binsPbPb[1];
+
+
+ for(Long_t iBin=1; iBin<= looppt;iBin++){
+ if(fBeamType==0)
+ {
+ nBinpp[0]=iBin;
+ backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
+ backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
+ }
+ if(fBeamType==1)
+ {
+ for(Long_t iiBin=1; iiBin<=binsPbPb[0];iiBin++){
+ nBinPbPb[0]=iiBin;
+ nBinPbPb[1]=iBin;
+ Double_t evtnormPbPb=0;
+ if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
+ backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
+ backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
+ }
+ }
+ }
+ //end of workaround for statistical errors
+ AliCFDataGrid *weightedNonHFEContainer;
+ if(fBeamType==0) weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,binspp);
+ else weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",2,binsPbPb);
weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
- CorrectStatErr(backgroundGrid);
+ delete[] nBinpp;
+ delete[] binspp;
+ delete[] nBinPbPb;
+ delete[] binsPbPb;
return backgroundGrid;
}
AliError("Data Container not available");
return NULL;
}
-
dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
}
-
AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
result->SetName("ParametrizedEfficiencyBefore");
THnSparse *h = result->GetGrid();
Int_t nbdimensions = h->GetNdimensions();
//printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
-
AliCFContainer *dataContainer = GetContainer(kDataContainer);
if(!dataContainer){
AliError("Data Container not available");
memset(coord, 0, sizeof(Int_t) * nbdimensions);
Double_t* points = new Double_t[nbdimensions];
-
ULong64_t nEntries = h->GetNbins();
for (ULong64_t i = 0; i < nEntries; ++i) {
AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
if(fDebugLevel > 0) {
-
+
TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
cEfficiencyParametrized->Divide(2,1);
cEfficiencyParametrized->cd(1);
//ratioefficiency->Divide(afterE);
//ratioefficiency->Draw();
+ if(fWriteToFile) cEfficiencyParametrized->SaveAs("efficiency.eps");
}
-
return resultt;
}
Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
for(Int_t binc = 0; binc < nbbin; binc++){
- TString titlee("V0Efficiency_centrality_bin_");
- titlee += binc;
- TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
- ccV0Efficiency->Divide(2,1);
- ccV0Efficiency->cd(1);
- gPad->SetLogy();
- cenaxisa->SetRange(binc+1,binc+1);
- cenaxisb->SetRange(binc+1,binc+1);
- cenaxisc->SetRange(binc+1,binc+1);
- TH1D *aftere = (TH1D *) sparseafter->Projection(1);
- TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
- CorrectFromTheWidth(aftere);
- CorrectFromTheWidth(beforee);
- aftere->SetStats(0);
- aftere->SetTitle((const char*)titlee);
- aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
- aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- aftere->SetMarkerStyle(25);
- aftere->SetMarkerColor(kBlack);
- aftere->SetLineColor(kBlack);
- beforee->SetStats(0);
- beforee->SetTitle((const char*)titlee);
- beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
- beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- beforee->SetMarkerStyle(24);
- beforee->SetMarkerColor(kBlue);
- beforee->SetLineColor(kBlue);
- aftere->Draw();
- beforee->Draw("same");
- TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
- lega->AddEntry(beforee,"Before correction","p");
- lega->AddEntry(aftere,"After correction","p");
- lega->Draw("same");
- cV0Efficiencye->cd(1);
- gPad->SetLogy();
- TH1D *afteree = (TH1D *) aftere->Clone();
- afteree->SetMarkerStyle(stylee[binc]);
- afteree->SetMarkerColor(colorr[binc]);
- if(binc==0) afteree->Draw();
- else afteree->Draw("same");
- legtotal->AddEntry(afteree,(const char*) titlee,"p");
- cV0Efficiencye->cd(2);
- gPad->SetLogy();
- TH1D *aftereeu = (TH1D *) aftere->Clone();
- aftereeu->SetMarkerStyle(stylee[binc]);
- aftereeu->SetMarkerColor(colorr[binc]);
- if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
- if(binc==0) aftereeu->Draw();
- else aftereeu->Draw("same");
- legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
- ccV0Efficiency->cd(2);
- TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
- efficiencyDDproj->SetTitle("");
- efficiencyDDproj->SetStats(0);
- efficiencyDDproj->SetMarkerStyle(25);
- efficiencyDDproj->Draw();
-
+ TString titlee("V0Efficiency_centrality_bin_");
+ titlee += binc;
+ TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
+ ccV0Efficiency->Divide(2,1);
+ ccV0Efficiency->cd(1);
+ gPad->SetLogy();
+ cenaxisa->SetRange(binc+1,binc+1);
+ cenaxisb->SetRange(binc+1,binc+1);
+ cenaxisc->SetRange(binc+1,binc+1);
+ TH1D *aftere = (TH1D *) sparseafter->Projection(1);
+ TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
+ CorrectFromTheWidth(aftere);
+ CorrectFromTheWidth(beforee);
+ aftere->SetStats(0);
+ aftere->SetTitle((const char*)titlee);
+ aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+ aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+ aftere->SetMarkerStyle(25);
+ aftere->SetMarkerColor(kBlack);
+ aftere->SetLineColor(kBlack);
+ beforee->SetStats(0);
+ beforee->SetTitle((const char*)titlee);
+ beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+ beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+ beforee->SetMarkerStyle(24);
+ beforee->SetMarkerColor(kBlue);
+ beforee->SetLineColor(kBlue);
+ aftere->Draw();
+ beforee->Draw("same");
+ TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
+ lega->AddEntry(beforee,"Before correction","p");
+ lega->AddEntry(aftere,"After correction","p");
+ lega->Draw("same");
+ cV0Efficiencye->cd(1);
+ gPad->SetLogy();
+ TH1D *afteree = (TH1D *) aftere->Clone();
+ afteree->SetMarkerStyle(stylee[binc]);
+ afteree->SetMarkerColor(colorr[binc]);
+ if(binc==0) afteree->Draw();
+ else afteree->Draw("same");
+ legtotal->AddEntry(afteree,(const char*) titlee,"p");
+ cV0Efficiencye->cd(2);
+ gPad->SetLogy();
+ TH1D *aftereeu = (TH1D *) aftere->Clone();
+ aftereeu->SetMarkerStyle(stylee[binc]);
+ aftereeu->SetMarkerColor(colorr[binc]);
+ if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
+ if(binc==0) aftereeu->Draw();
+ else aftereeu->Draw("same");
+ legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
+ ccV0Efficiency->cd(2);
+ TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
+ efficiencyDDproj->SetTitle("");
+ efficiencyDDproj->SetStats(0);
+ efficiencyDDproj->SetMarkerStyle(25);
+ efficiencyDDproj->Draw();
}
cV0Efficiencye->cd(1);
cenaxisa->SetRange(0,nbbin);
cenaxisb->SetRange(0,nbbin);
cenaxisc->SetRange(0,nbbin);
-
+
+ if(fWriteToFile) cV0Efficiencye->SaveAs("V0efficiency.eps");
}
}
// Efficiency
AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
-
- // Consider parameterized IP cut efficiency
- if(!fInclusiveSpectrum){
- Int_t* bins=new Int_t[1];
- bins[0]=44;
- AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
- beffContainer->SetGrid(GetBeautyIPEff());
- efficiencyD->Multiply(beffContainer,1);
+ if(!fBeauty2ndMethod)
+ {
+ // Consider parameterized IP cut efficiency
+ if(!fInclusiveSpectrum){
+ Int_t* bins=new Int_t[1];
+ bins[0]=35;
+
+ AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
+ beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
+ efficiencyD->Multiply(beffContainer,1);
+ }
}
// Unfold
AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
- //unfolding.SetUseCorrelatedErrors();
if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
if(fSetSmoothing) unfolding.UseSmoothing();
ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
ratioresidual->SetStats(0);
ratioresidual->Draw();
-
+
+ if(fWriteToFile) cresidual->SaveAs("Unfolding.eps");
}
return listofresults;
AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
- // Consider parameterized IP cut efficiency
- if(!fInclusiveSpectrum){
- Int_t* bins=new Int_t[1];
- bins[0]=44;
-
- AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
- beffContainer->SetGrid(GetBeautyIPEff());
- efficiencyD->Multiply(beffContainer,1);
+
+ if(!fBeauty2ndMethod)
+ {
+ // Consider parameterized IP cut efficiency
+ if(!fInclusiveSpectrum){
+ Int_t* bins=new Int_t[1];
+ bins[0]=35;
+
+ AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
+ beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
+ efficiencyD->Multiply(beffContainer,1);
+ }
}
// Data in the right format
//Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
//Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
- TString titlee("Efficiency_centrality_bin_");
- titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
- titlee += "_";
- titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
- TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
- cefficiencye->Divide(2,1);
- cefficiencye->cd(1);
- gPad->SetLogy();
- cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
- cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
- TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
- TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
- CorrectFromTheWidth(afterefficiency);
- CorrectFromTheWidth(beforeefficiency);
- afterefficiency->SetStats(0);
- afterefficiency->SetTitle((const char*)titlee);
- afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
- afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- afterefficiency->SetMarkerStyle(25);
- afterefficiency->SetMarkerColor(kBlack);
- afterefficiency->SetLineColor(kBlack);
- beforeefficiency->SetStats(0);
- beforeefficiency->SetTitle((const char*)titlee);
- beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
- beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- beforeefficiency->SetMarkerStyle(24);
- beforeefficiency->SetMarkerColor(kBlue);
- beforeefficiency->SetLineColor(kBlue);
- afterefficiency->Draw();
- beforeefficiency->Draw("same");
- TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
- lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
- lega->AddEntry(afterefficiency,"After efficiency correction","p");
- lega->Draw("same");
- cefficiencye->cd(2);
- cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
- TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
- efficiencyDDproj->SetTitle("");
- efficiencyDDproj->SetStats(0);
- efficiencyDDproj->SetMarkerStyle(25);
- efficiencyDDproj->SetMarkerColor(2);
- efficiencyDDproj->Draw();
- cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
- TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
- efficiencyDDproja->SetTitle("");
- efficiencyDDproja->SetStats(0);
- efficiencyDDproja->SetMarkerStyle(26);
- efficiencyDDproja->SetMarkerColor(4);
- efficiencyDDproja->Draw("same");
-
+ TString titlee("Efficiency_centrality_bin_");
+ titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
+ titlee += "_";
+ titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
+ TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
+ cefficiencye->Divide(2,1);
+ cefficiencye->cd(1);
+ gPad->SetLogy();
+ cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
+ cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
+ TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
+ TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
+ CorrectFromTheWidth(afterefficiency);
+ CorrectFromTheWidth(beforeefficiency);
+ afterefficiency->SetStats(0);
+ afterefficiency->SetTitle((const char*)titlee);
+ afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+ afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+ afterefficiency->SetMarkerStyle(25);
+ afterefficiency->SetMarkerColor(kBlack);
+ afterefficiency->SetLineColor(kBlack);
+ beforeefficiency->SetStats(0);
+ beforeefficiency->SetTitle((const char*)titlee);
+ beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+ beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+ beforeefficiency->SetMarkerStyle(24);
+ beforeefficiency->SetMarkerColor(kBlue);
+ beforeefficiency->SetLineColor(kBlue);
+ afterefficiency->Draw();
+ beforeefficiency->Draw("same");
+ TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
+ lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
+ lega->AddEntry(afterefficiency,"After efficiency correction","p");
+ lega->Draw("same");
+ cefficiencye->cd(2);
+ cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
+ TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
+ efficiencyDDproj->SetTitle("");
+ efficiencyDDproj->SetStats(0);
+ efficiencyDDproj->SetMarkerStyle(25);
+ efficiencyDDproj->SetMarkerColor(2);
+ efficiencyDDproj->Draw();
+ cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
+ TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
+ efficiencyDDproja->SetTitle("");
+ efficiencyDDproja->SetStats(0);
+ efficiencyDDproja->SetMarkerStyle(26);
+ efficiencyDDproja->SetMarkerColor(4);
+ efficiencyDDproja->Draw("same");
}
}
-
+ if(fWriteToFile) cefficiency->SaveAs("efficiencyCorrected.eps");
}
return result;
// Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
// Give the final pt spectrum to be compared
//
-
+
if(fNEvents[i] > 0) {
Int_t ptpr = 0;
TH1D* projection = (TH1D *) spectrum->Project(ptpr);
CorrectFromTheWidth(projection);
TGraphErrors *graphError = NormalizeTH1(projection,i);
+
return graphError;
}
// Give the final pt spectrum to be compared
//
Double_t chargecoefficient = 0.5;
- if(fChargeChoosen >= 0) chargecoefficient = 1.0;
+ if(fChargeChoosen != 0) chargecoefficient = 1.0;
- Double_t etarange = fEtaSelected ? fEtaRange[1] - fEtaRange[0] : 1.6;
+ Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
printf("Normalizing Eta Range %f\n", etarange);
if(fNEvents[i] > 0) {
dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
n = input->GetBinContent(ibin);
dN = input->GetBinError(ibin);
-
// New point
nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
errdN = 1./(2. * TMath::Pi() * p);
spectrumNormalized->SetMarkerStyle(22);
spectrumNormalized->SetMarkerColor(kBlue);
spectrumNormalized->SetLineColor(kBlue);
-
+
return spectrumNormalized;
}
// Give the final pt spectrum to be compared
//
Double_t chargecoefficient = 0.5;
- if(fChargeChoosen >= 0) chargecoefficient = 1.0;
+ if(fChargeChoosen != kAllCharge) chargecoefficient = 1.0;
- Double_t etarange = fEtaSelected ? fEtaRange[1] - fEtaRange[0] : 1.6;
+ Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
printf("Normalizing Eta Range %f\n", etarange);
if(normalization > 0) {
dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
n = input->GetBinContent(ibin);
dN = input->GetBinError(ibin);
-
// New point
nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
errdN = 1./(2. * TMath::Pi() * p);
//
// Set the container for a given step to the
//
- if(!fCFContainers) fCFContainers = new TList;
+ if(!fCFContainers) fCFContainers = new TObjArray(kDataContainerV0+1);
fCFContainers->AddAt(cont, type);
}
return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
}
//____________________________________________________________
-AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Int_t positivenegative) {
+AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh) {
//
// Slice bin for a given source of electron
// nDim is the number of dimension the corrections are done
// dimensions are the definition of the dimensions
// source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
// positivenegative if we want to keep positive (1) or negative (0) or both (-1)
+ // centrality (-1 means we do not cut on centrality)
//
Double_t *varMin = new Double_t[container->GetNVar()],
// source
if(ivar == 4){
if((source>= 0) && (source<container->GetNBins(ivar))) {
- varMin[ivar] = binLimits[source];
- varMax[ivar] = binLimits[source];
+ varMin[ivar] = binLimits[source];
+ varMax[ivar] = binLimits[source];
}
}
// charge
if(ivar == 3) {
- if((positivenegative>= 0) && (positivenegative<container->GetNBins(ivar))) {
- varMin[ivar] = binLimits[positivenegative];
- varMax[ivar] = binLimits[positivenegative];
- }
+ if(charge != kAllCharge) varMin[ivar] = varMax[ivar] = charge;
}
// eta
if(ivar == 1){
- for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++) AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic)));
+ for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++)
+ AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic)));
if(fEtaSelected){
varMin[ivar] = fEtaRange[0];
varMax[ivar] = fEtaRange[1];
}
}
-
+ if(fEtaSelected){
+ fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
+ fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
+ AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0]));
+ }
+ if(ivar == 5){
+ if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar))) {
+ varMin[ivar] = binLimits[centralitylow];
+ varMax[ivar] = binLimits[centralityhigh];
+
+ TAxis *axistest = container->GetAxis(5,0);
+ printf("GetSlicedContainer: Number of bin in centrality direction %d\n",axistest->GetNbins());
+ printf("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]);
+ Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
+ Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
+ printf("GetSlicedContainer: Low centrality %f and high centrality %f\n",lowcentrality,highcentrality);
+
+ }
+
+ }
+
delete[] binLimits;
}
}
//_________________________________________________________________________
-THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions) const {
+THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions, Int_t centralitylow, Int_t centralityhigh) const {
//
// Slice correlation
//
AliError("Problem in the dimensions");
return NULL;
}
+
+ // Cut in centrality is centrality > -1
+ if((centralitylow >=0) && (centralityhigh >=0)) {
+
+ TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
+ TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
+
+ Int_t bins0 = axiscentrality0->GetNbins();
+ Int_t bins1 = axiscentrality1->GetNbins();
+
+ printf("Number of centrality bins: %d and %d\n",bins0,bins1);
+ if(bins0 != bins1) {
+ AliError("Problem in the dimensions");
+ return NULL;
+ }
+
+ if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0)) {
+ axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
+ axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
+
+ Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
+ Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
+ Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
+ Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
+ printf("GetCorrelation0: Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0);
+ printf("GetCorrelation1: Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1);
+
+ TCanvas * ctestcorrelation = new TCanvas("testcorrelationprojection","testcorrelationprojection",1000,700);
+ ctestcorrelation->cd(1);
+ TH2D* projection = (TH2D *) correlationmatrix->Projection(5,5+((Int_t)(ndimensions/2.)));
+ projection->Draw("colz");
+
+ }
+
+ }
+
+
Int_t ndimensionsContainer = (Int_t) ndimensions/2;
//printf("Number of dimension %d container\n",ndimensionsContainer);
// Measured D->e based weighting factors
//
- const Int_t nDim=1;
- Int_t nBin[nDim] = {44};
- const Double_t kPtbound[2] = {0.1, 20.};
+ const Int_t nDimpp=1;
+ Int_t nBinpp[nDimpp] = {35};
+ Double_t ptbinning1[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
+ const Int_t nDimPbPb=2;
+ Int_t nBinPbPb[nDimPbPb] = {11,35};
+ Double_t kCentralityRange[12] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
+ Int_t loopcentr=1;
+ Int_t looppt=nBinpp[0];
+ if(fBeamType==0)
+ {
+ fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
+ fWeightCharm->SetBinEdges(0, ptbinning1);
+ }
+ if(fBeamType==1)
+ {
+ fWeightCharm = new THnSparseF("weightHisto", "weighting factor; centrality bin; pt[GeV/c]", nDimPbPb, nBinPbPb);
+ fWeightCharm->SetBinEdges(1, ptbinning1);
+ fWeightCharm->SetBinEdges(0, kCentralityRange);
+ loopcentr=nBinPbPb[0];
+ }
+
+ // Weighting factor for pp
+ Double_t weight[35]={0.859260, 0.872552, 0.847475, 0.823631, 0.839386, 0.874024, 0.916755, 0.942801, 0.965856, 0.933905, 0.933414, 0.931936, 0.847826, 0.810902, 0.796608, 0.727002, 0.659227, 0.583610, 0.549956, 0.512633, 0.472254, 0.412364, 0.353191, 0.319145, 0.305412, 0.290334, 0.269863, 0.254646, 0.230245, 0.200859, 0.275953, 0.276271, 0.227332, 0.197004, 0.474385};
+
+ // Weighting factor for PbPb (0-20%)
+ //Double_t weight[35]={0.641897, 0.640472, 0.615228, 0.650469, 0.737762, 0.847867, 1.009317, 1.158594, 1.307482, 1.476973, 1.551131, 1.677131, 1.785478, 1.888933, 2.017957, 2.074757, 1.926700, 1.869495, 1.546558, 1.222873, 1.160313, 0.903375, 0.799642, 0.706244, 0.705449, 0.599947, 0.719570, 0.499422, 0.703978, 0.477452, 0.325057, 0.093391, 0.096675, 0.000000, 0.000000};
- Double_t* binEdges[nDim];
- binEdges[0] = AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
+ // Weighting factor for PbPb (40-80%)
+ //Double_t weight[35]={0.181953, 0.173248, 0.166799, 0.182558, 0.206581, 0.236955, 0.279390, 0.329129, 0.365260, 0.423059, 0.452057, 0.482726, 0.462627, 0.537770, 0.584663, 0.579452, 0.587194, 0.499498, 0.443299, 0.398596, 0.376695, 0.322331, 0.260890, 0.374834, 0.249114, 0.310330, 0.287326, 0.243174, 0.758945, 0.138867, 0.170576, 0.107797, 0.011390, 0.000000, 0.000000};
- fWeightCharm = new THnSparseF("weightHisto", "weiting factor; pt[GeV/c]", nDim, nBin);
- for(Int_t idim = 0; idim < nDim; idim++)
- fWeightCharm->SetBinEdges(idim, binEdges[idim]);
- Double_t weight[44]={
-1.11865,1.11444,1.13395,1.15751,1.13412,1.14446,1.15372,1.09458,1.13446,1.11707,1.1352,1.10979,1.08887,1.11164,1.10772,1.10727,1.07027,1.07016,1.04826,1.03248,1.0093,0.973534,0.932574,0.869838,0.799316,0.734015,0.643001,0.584319,0.527147,0.495661,0.470002,0.441129,0.431156,0.417242,0.39246,0.379611,0.390845,0.368857,0.34959,0.387186,0.376364,0.384437,0.349585,0.383225};
//points
Double_t pt[1];
- for(int i=0; i<nBin[0]; i++){
- pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
- fWeightCharm->Fill(pt,weight[i]);
+ Double_t contents[2];
+
+ for(int icentr=0; icentr<loopcentr; icentr++)
+ {
+ for(int i=0; i<looppt; i++){
+ pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
+ if(fBeamType==1)
+ {
+ contents[0]=icentr;
+ contents[1]=pt[0];
+ }
+ if(fBeamType==0)
+ {
+ contents[0]=pt[0];
+ }
+
+ fWeightCharm->Fill(contents,weight[i]);
+ }
}
- Int_t* ibins[nDim];
- for(Int_t ivar = 0; ivar < nDim; ivar++)
- ibins[ivar] = new Int_t[nBin[ivar] + 1];
- fWeightCharm->SetBinError(ibins[0],0);
+
+ Int_t nDimSparse = fWeightCharm->GetNdimensions();
+ Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
+ Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
+
+ for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+ binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
+ nBins *= binsvar[iVar];
+ }
+
+ Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+ // loop that sets 0 error in each bin
+ for (Long_t iBin=0; iBin<nBins; iBin++) {
+ Long_t bintmp = iBin ;
+ for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+ binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+ bintmp /= binsvar[iVar] ;
+ }
+ fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
+ }
+
+ delete[] binsvar;
+ delete[] binfill;
return fWeightCharm;
}
//____________________________________________________________________________
-THnSparse* AliHFEspectrum::GetBeautyIPEff(){
+void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
+
+ // TOF PID efficiencies
+ Int_t ptpr;
+ if(fBeamType==0) ptpr=0;
+ if(fBeamType==1) ptpr=1;
+
+ Int_t loopcentr=1;
+ const Int_t nCentralitybinning=11; //number of centrality bins
+ if(fBeamType==1)
+ {
+ loopcentr=nCentralitybinning;
+ }
+
+ TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
+ TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
+ TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.3,10.0);
+
+ TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
+ cefficiencyParamtof->cd();
+
+ AliCFContainer *mccontainermcD = 0x0;
+ AliCFContainer *mccontaineresdD = 0x0;
+ TH1D* efficiencysigTOFPIDD[nCentralitybinning];
+ TH1D* efficiencyTOFPIDD[nCentralitybinning];
+ TH1D* efficiencysigesdTOFPIDD[nCentralitybinning];
+ TH1D* efficiencyesdTOFPIDD[nCentralitybinning];
+ Int_t source = -1; //get parameterized TOF PID efficiencies
+
+ for(int icentr=0; icentr<loopcentr; icentr++) {
+ // signal sample
+ if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
+ AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
+ efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
+
+ // mb sample for double check
+ if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
+ AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
+ efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
+
+ // mb sample with reconstructed variables
+ if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
+ AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
+ efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
+
+ // mb sample with reconstructed variables
+ if(fBeamType==0) mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
+ AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
+ efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
+
+ //fill histo
+ efficiencysigTOFPIDD[icentr] = (TH1D *) efficiencymcsigParamTOFPID->Project(ptpr);
+ efficiencyTOFPIDD[icentr] = (TH1D *) efficiencymcParamTOFPID->Project(ptpr);
+ efficiencysigesdTOFPIDD[icentr] = (TH1D *) efficiencysigesdParamTOFPID->Project(ptpr);
+ efficiencyesdTOFPIDD[icentr] = (TH1D *) efficiencyesdParamTOFPID->Project(ptpr);
+ efficiencysigTOFPIDD[icentr]->SetName(Form("efficiencysigTOFPIDD%d",icentr));
+ efficiencyTOFPIDD[icentr]->SetName(Form("efficiencyTOFPIDD%d",icentr));
+ efficiencysigesdTOFPIDD[icentr]->SetName(Form("efficiencysigesdTOFPIDD%d",icentr));
+ efficiencyesdTOFPIDD[icentr]->SetName(Form("efficiencyesdTOFPIDD%d",icentr));
+
+ //fit (mc enhenced sample)
+ fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+ efficiencysigTOFPIDD[icentr]->Fit(fittofpid,"R");
+ efficiencysigTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
+ fEfficiencyTOFPIDD[icentr] = efficiencysigTOFPIDD[icentr]->GetFunction("fittofpid");
+
+ //fit (esd enhenced sample)
+ efficiencysigesdTOFPIDD[icentr]->Fit(fittofpid,"R");
+ efficiencysigesdTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
+ fEfficiencyesdTOFPIDD[icentr] = efficiencysigesdTOFPIDD[icentr]->GetFunction("fittofpid");
+
+ }
+
+ // draw (for PbPb, only 1st bin)
+ //sig mc
+ efficiencysigTOFPIDD[0]->SetTitle("");
+ efficiencysigTOFPIDD[0]->SetStats(0);
+ efficiencysigTOFPIDD[0]->SetMarkerStyle(25);
+ efficiencysigTOFPIDD[0]->SetMarkerColor(2);
+ efficiencysigTOFPIDD[0]->SetLineColor(2);
+ efficiencysigTOFPIDD[0]->Draw();
+
+ //mb mc
+ efficiencyTOFPIDD[0]->SetTitle("");
+ efficiencyTOFPIDD[0]->SetStats(0);
+ efficiencyTOFPIDD[0]->SetMarkerStyle(24);
+ efficiencyTOFPIDD[0]->SetMarkerColor(4);
+ efficiencyTOFPIDD[0]->SetLineColor(4);
+ efficiencyTOFPIDD[0]->Draw("same");
+
+ //sig esd
+ efficiencysigesdTOFPIDD[0]->SetTitle("");
+ efficiencysigesdTOFPIDD[0]->SetStats(0);
+ efficiencysigesdTOFPIDD[0]->SetMarkerStyle(25);
+ efficiencysigesdTOFPIDD[0]->SetMarkerColor(3);
+ efficiencysigesdTOFPIDD[0]->SetLineColor(3);
+ efficiencysigesdTOFPIDD[0]->Draw("same");
+
+ //mb esd
+ efficiencyesdTOFPIDD[0]->SetTitle("");
+ efficiencyesdTOFPIDD[0]->SetStats(0);
+ efficiencyesdTOFPIDD[0]->SetMarkerStyle(25);
+ efficiencyesdTOFPIDD[0]->SetMarkerColor(1);
+ efficiencyesdTOFPIDD[0]->SetLineColor(1);
+ efficiencyesdTOFPIDD[0]->Draw("same");
+
+ //signal mc fit
+ if(fEfficiencyTOFPIDD[0]){
+ fEfficiencyTOFPIDD[0]->SetLineColor(2);
+ fEfficiencyTOFPIDD[0]->Draw("same");
+ }
+ //mb esd fit
+ if(fEfficiencyesdTOFPIDD[0]){
+ fEfficiencyesdTOFPIDD[0]->SetLineColor(3);
+ fEfficiencyesdTOFPIDD[0]->Draw("same");
+ }
+
+ TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
+ legtofeff->AddEntry(efficiencysigTOFPIDD[0],"TOF PID Step Efficiency","");
+ legtofeff->AddEntry(efficiencysigTOFPIDD[0],"vs MC p_{t} for enhenced samples","p");
+ legtofeff->AddEntry(efficiencyTOFPIDD[0],"vs MC p_{t} for mb samples","p");
+ legtofeff->AddEntry(efficiencysigesdTOFPIDD[0],"vs esd p_{t} for enhenced samples","p");
+ legtofeff->AddEntry(efficiencyesdTOFPIDD[0],"vs esd p_{t} for mb samples","p");
+ legtofeff->Draw("same");
+
+
+ TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
+ cefficiencyParamIP->cd();
+ gStyle->SetOptStat(0);
+
+ // IP cut efficiencies
+ for(int icentr=0; icentr<loopcentr; icentr++) {
+
+ AliCFContainer *charmCombined = 0x0;
+ AliCFContainer *beautyCombined = 0x0;
+ AliCFContainer *beautyCombinedesd = 0x0;
+
+ printf("centrality printing %i \n",icentr);
+
+ source = 0; //charm enhenced
+ if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+ AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
+ efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+ charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");
+
+ source = 1; //beauty enhenced
+ if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+ AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
+ efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+ beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined");
+
+ if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+ AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
+ efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+ beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
+
+ source = 0; //charm mb
+ if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+ AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
+ efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+ charmCombined->Add(mccontainermcD);
+ AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined);
+ efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
+
+ source = 1; //beauty mb
+ if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+ AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
+ efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+ beautyCombined->Add(mccontainermcD);
+ AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined);
+ efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
+
+ if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+ AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
+ efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+ beautyCombinedesd->Add(mccontaineresdD);
+ AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
+ efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
+
+ source = 2; //conversion mb
+ if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+ AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
+ efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+ source = 3; //non HFE except for the conversion mb
+ if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
+ else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+ AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
+ efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+ if(fIPEffCombinedSamples){
+ fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmCombined->Project(ptpr); //signal enhenced + mb
+ fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautyCombined->Project(ptpr); //signal enhenced + mb
+ fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautyCombinedesd->Project(ptpr); //signal enhenced + mb
+ }
+ else{
+ fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr); //signal enhenced only
+ fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr); //signal enhenced only
+ fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautySigesd->Project(ptpr); //signal enhenced only
+ }
+ fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr); //mb only
+ fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr); //mb only
+ fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr); //mb only
+ fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr); //mb only
+
+ }
+
+ if(fBeamType==0){
+ AliCFEffGrid *nonHFEEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
+ fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
+
+ AliCFEffGrid *conversionEffGrid = (AliCFEffGrid*) GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
+ fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
+ }
+
+ for(int icentr=0; icentr<loopcentr; icentr++) {
+ fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+ fipfit->SetLineColor(2);
+ fEfficiencyBeautySigD[icentr]->Fit(fipfit,"R");
+ fEfficiencyBeautySigD[icentr]->GetYaxis()->SetTitle("Efficiency");
+ if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit"); //why do we need this line?
+ else fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[icentr]->GetFunction("fipfit");
+
+ fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+ fipfit->SetLineColor(6);
+ fEfficiencyBeautySigesdD[icentr]->Fit(fipfit,"R");
+ fEfficiencyBeautySigesdD[icentr]->GetYaxis()->SetTitle("Efficiency");
+ if(fBeauty2ndMethod)fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[0]->GetFunction("fipfit"); //why do we need this line?
+ else fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[icentr]->GetFunction("fipfit");
+
+ fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+ fipfit->SetLineColor(1);
+ fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R");
+ fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency");
+ fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit");
+
+ if(fIPParameterizedEff){
+ fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+ fipfitnonhfe->SetLineColor(3);
+ fConversionEff[icentr]->Fit(fipfitnonhfe,"R");
+ fConversionEff[icentr]->GetYaxis()->SetTitle("Efficiency");
+ fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfitnonhfe");
+
+ fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+ fipfitnonhfe->SetLineColor(4);
+ fNonHFEEff[icentr]->Fit(fipfitnonhfe,"R");
+ fNonHFEEff[icentr]->GetYaxis()->SetTitle("Efficiency");
+ fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfitnonhfe");
+ }
+ }
+
+ // draw (for PbPb, only 1st bin)
+ fEfficiencyCharmSigD[0]->SetMarkerStyle(21);
+ fEfficiencyCharmSigD[0]->SetMarkerColor(1);
+ fEfficiencyCharmSigD[0]->SetLineColor(1);
+ fEfficiencyBeautySigD[0]->SetMarkerStyle(21);
+ fEfficiencyBeautySigD[0]->SetMarkerColor(2);
+ fEfficiencyBeautySigD[0]->SetLineColor(2);
+ fEfficiencyBeautySigesdD[0]->SetStats(0);
+ fEfficiencyBeautySigesdD[0]->SetMarkerStyle(21);
+ fEfficiencyBeautySigesdD[0]->SetMarkerColor(6);
+ fEfficiencyBeautySigesdD[0]->SetLineColor(6);
+ fCharmEff[0]->SetMarkerStyle(24);
+ fCharmEff[0]->SetMarkerColor(1);
+ fCharmEff[0]->SetLineColor(1);
+ fBeautyEff[0]->SetMarkerStyle(24);
+ fBeautyEff[0]->SetMarkerColor(2);
+ fBeautyEff[0]->SetLineColor(2);
+ fConversionEff[0]->SetMarkerStyle(24);
+ fConversionEff[0]->SetMarkerColor(3);
+ fConversionEff[0]->SetLineColor(3);
+ fNonHFEEff[0]->SetMarkerStyle(24);
+ fNonHFEEff[0]->SetMarkerColor(4);
+ fNonHFEEff[0]->SetLineColor(4);
+
+ fEfficiencyCharmSigD[0]->Draw();
+ fEfficiencyCharmSigD[0]->GetXaxis()->SetRangeUser(0.0,7.9);
+ fEfficiencyCharmSigD[0]->GetYaxis()->SetRangeUser(0.0,0.5);
+
+ fEfficiencyBeautySigD[0]->Draw("same");
+ fEfficiencyBeautySigesdD[0]->Draw("same");
+ //fCharmEff[0]->Draw("same");
+ //fBeautyEff[0]->Draw("same");
+
+ if(fBeamType==0){
+ fConversionEffbgc->SetMarkerStyle(25);
+ fConversionEffbgc->SetMarkerColor(3);
+ fConversionEffbgc->SetLineColor(3);
+ fNonHFEEffbgc->SetMarkerStyle(25);
+ fNonHFEEffbgc->SetMarkerColor(4);
+ fNonHFEEffbgc->SetLineColor(4);
+ fConversionEffbgc->Draw("same");
+ fNonHFEEffbgc->Draw("same");
+ }
+ else{
+ fConversionEff[0]->Draw("same");
+ fNonHFEEff[0]->Draw("same");
+ }
+ if(fEfficiencyIPBeautyD[0])
+ fEfficiencyIPBeautyD[0]->Draw("same");
+ if(fEfficiencyIPBeautyesdD[0])
+ fEfficiencyIPBeautyesdD[0]->Draw("same");
+ if( fEfficiencyIPCharmD[0])
+ fEfficiencyIPCharmD[0]->Draw("same");
+ if(fIPParameterizedEff){
+ fEfficiencyIPConversionD[0]->Draw("same");
+ fEfficiencyIPNonhfeD[0]->Draw("same");
+ }
+ TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
+ legipeff->AddEntry(fEfficiencyBeautySigD[0],"IP Step Efficiency","");
+ legipeff->AddEntry(fEfficiencyBeautySigD[0],"beauty e","p");
+ legipeff->AddEntry(fEfficiencyBeautySigesdD[0],"beauty e(esd pt)","p");
+ legipeff->AddEntry(fEfficiencyCharmSigD[0],"charm e","p");
+ legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
+ legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
+ //legipeff->AddEntry(fConversionEff[0],"conversion e","p");
+ //legipeff->AddEntry(fNonHFEEff[0],"Dalitz e","p");
+ legipeff->Draw("same");
+ gPad->SetGrid();
+ //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
+}
+
+//____________________________________________________________________________
+THnSparse* AliHFEspectrum::GetBeautyIPEff(Bool_t isMCpt){
//
// Return beauty electron IP cut efficiency
//
- const Int_t nDim=1;
- Int_t nBin[nDim] = {44};
- const Double_t kPtbound[2] = {0.1, 20.};
+ const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
+ const Int_t nCentralitybinning=11;//number of centrality bins
+ Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
+ Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
+ Int_t ptpr = 0;
+ Int_t nDim=1; //dimensions of the efficiency weighting grid
+ if(fBeamType==1)
+ {
+ ptpr=1;
+ nDim=2; //dimensions of the efficiency weighting grid
+ }
+ Int_t nBin[1] = {nPtbinning1};
+ Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
- Double_t* binEdges[nDim];
- binEdges[0] = AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
- THnSparseF *ipcut = new THnSparseF("beff", "b eff; pt[GeV/c]", nDim, nBin);
+ THnSparseF *ipcut;
+ if(fBeamType==0) ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
+ else ipcut = new THnSparseF("beff", "b IP efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
+
for(Int_t idim = 0; idim < nDim; idim++)
- ipcut->SetBinEdges(idim, binEdges[idim]);
+ {
+ if(nDim==1) ipcut->SetBinEdges(idim, kPtRange);
+ if(nDim==2)
+ {
+ ipcut->SetBinEdges(0, kCentralityRange);
+ ipcut->SetBinEdges(1, kPtRange);
+ }
+ }
Double_t pt[1];
- Double_t weight;
- for(int i=0; i<nBin[0]; i++){
- pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
- weight=0.612*(0.385+0.111*log(pt[0])-0.00435*log(pt[0])*log(pt[0]))*tanh(5.42*pt[0]-1.22); // for 3 sigma cut
- //weight=0.786*(0.493+0.0283*log(pt[0])-0.0190*log(pt[0])*log(pt[0]))*tanh(1.84*pt[0]-0.00368); // for 2 sigma cut
- ipcut->Fill(pt,weight);
- }
- Int_t* ibins[nDim];
- for(Int_t ivar = 0; ivar < nDim; ivar++)
- ibins[ivar] = new Int_t[nBin[ivar] + 1];
- ipcut->SetBinError(ibins[0],0);
+ Double_t weight[nCentralitybinning];
+ Double_t weightErr[nCentralitybinning];
+ Double_t contents[2];
+
+ for(Int_t a=0;a<11;a++)
+ {
+ weight[a] = 1.0;
+ weightErr[a] = 1.0;
+ }
+
+
+ Int_t looppt=nBin[0];
+ Int_t loopcentr=1;
+ Int_t ibin[2];
+ if(fBeamType==1)
+ {
+ loopcentr=nBinPbPb[0];
+ }
+
+
+ for(int icentr=0; icentr<loopcentr; icentr++)
+ {
+ for(int i=0; i<looppt; i++)
+ {
+ pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
+ if(isMCpt){
+ if(fEfficiencyIPBeautyD[icentr]){
+ weight[icentr]=fEfficiencyIPBeautyD[icentr]->Eval(pt[0]);
+ weightErr[icentr] = 0;
+ }
+ else{
+ printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
+ weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1);
+ weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
+ }
+ }
+ else{
+ if(fEfficiencyIPBeautyesdD[icentr]){
+ weight[icentr]=fEfficiencyIPBeautyesdD[icentr]->Eval(pt[0]);
+ weightErr[icentr] = 0;
+ }
+ else{
+ printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
+ weight[icentr] = fEfficiencyBeautySigesdD[icentr]->GetBinContent(i+1);
+ weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
+ }
+ }
+
+ if(fBeamType==1){
+ contents[0]=icentr;
+ contents[1]=pt[0];
+ ibin[0]=icentr;
+ ibin[1]=i+1;
+ }
+ if(fBeamType==0){
+ contents[0]=pt[0];
+ ibin[0]=i+1;
+ }
+ ipcut->Fill(contents,weight[icentr]);
+ ipcut->SetBinError(ibin,weightErr[icentr]);
+ }
+ }
+
+ Int_t nDimSparse = ipcut->GetNdimensions();
+ Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
+ Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
+
+ for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+ binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
+ nBins *= binsvar[iVar];
+ }
+
+ Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+ // loop that sets 0 error in each bin
+ for (Long_t iBin=0; iBin<nBins; iBin++) {
+ Long_t bintmp = iBin ;
+ for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+ binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+ bintmp /= binsvar[iVar] ;
+ }
+ //ipcut->SetBinError(binfill,0.); // put 0 everywhere
+ }
+
+ delete[] binsvar;
+ delete[] binfill;
return ipcut;
}
//____________________________________________________________________________
-THnSparse* AliHFEspectrum::GetCharmEff(){
+THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
//
- // Return charm electron IP cut efficiency
+ // Return PID x IP cut efficiency
//
+ const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
+ const Int_t nCentralitybinning=11;//number of centrality bins
+ Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
+ Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
+ Int_t ptpr = 0;
+ Int_t nDim=1; //dimensions of the efficiency weighting grid
+ if(fBeamType==1)
+ {
+ ptpr=1;
+ nDim=2; //dimensions of the efficiency weighting grid
+ }
+ Int_t nBin[1] = {nPtbinning1};
+ Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
- const Int_t nDim=1;
- Int_t nBin[nDim] = {44};
- const Double_t kPtbound[2] = {0.1, 20.};
+ THnSparseF *pideff;
+ if(fBeamType==0) pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
+ else pideff = new THnSparseF("pideff", "PID efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
+ for(Int_t idim = 0; idim < nDim; idim++)
+ {
- Double_t* binEdges[nDim];
- binEdges[0] = AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
+ if(nDim==1) pideff->SetBinEdges(idim, kPtRange);
+ if(nDim==2)
+ {
+ pideff->SetBinEdges(0, kCentralityRange);
+ pideff->SetBinEdges(1, kPtRange);
+ }
+ }
- THnSparseF *ipcut = new THnSparseF("ceff", "c eff; pt[GeV/c]", nDim, nBin);
- for(Int_t idim = 0; idim < nDim; idim++)
- ipcut->SetBinEdges(idim, binEdges[idim]);
Double_t pt[1];
- Double_t weight;
- for(int i=0; i<nBin[0]; i++){
- pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
- weight=0.299*(0.307+0.176*log(pt[0])+0.0176*log(pt[0])*log(pt[0]))*tanh(96.3*pt[0]-19.7); // for 3 sigma cut
- //weight=0.477*(0.3757+0.184*log(pt[0])-0.0013*log(pt[0])*log(pt[0]))*tanh(436.013*pt[0]-96.2137); // for 2 sigma cut
- ipcut->Fill(pt,weight);
- }
- Int_t* ibins[nDim];
- for(Int_t ivar = 0; ivar < nDim; ivar++)
- ibins[ivar] = new Int_t[nBin[ivar] + 1];
- ipcut->SetBinError(ibins[0],0);
+ Double_t weight[nCentralitybinning];
+ Double_t weightErr[nCentralitybinning];
+ Double_t contents[2];
+
+ for(Int_t a=0;a<nCentralitybinning;a++)
+ {
+ weight[a] = 1.0;
+ weightErr[a] = 1.0;
+ }
+
+ Int_t looppt=nBin[0];
+ Int_t loopcentr=1;
+ Int_t ibin[2];
+ if(fBeamType==1)
+ {
+ loopcentr=nBinPbPb[0];
+ }
+
+ for(int icentr=0; icentr<loopcentr; icentr++)
+ {
+ Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
+ for(int i=0; i<looppt; i++)
+ {
+ pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
+
+ Double_t tofpideff = 0.;
+ Double_t tofpideffesd = 0.;
+ if(fEfficiencyTOFPIDD[icentr])
+ tofpideff = fEfficiencyTOFPIDD[icentr]->Eval(pt[0]);
+ else{
+ printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
+ }
+ if(fEfficiencyesdTOFPIDD[icentr])
+ tofpideffesd = fEfficiencyesdTOFPIDD[icentr]->Eval(pt[0]);
+ else{
+ printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
+ }
+
+ //tof pid eff x tpc pid eff x ip cut eff
+ if(fIPParameterizedEff){
+ if(source==0) {
+ if(fEfficiencyIPCharmD[icentr]){
+ weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
+ weightErr[icentr] = 0;
+ }
+ else{
+ printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
+ weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
+ weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
+ }
+ }
+ else if(source==2) {
+ if(fEfficiencyIPConversionD[icentr]){
+ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]);
+ weightErr[icentr] = 0;
+ }
+ else{
+ printf("Fit failed on conversion IP cut efficiency for centrality %d\n",icentr);
+ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1);
+ weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
+ }
+ }
+ else if(source==3) {
+ if(fEfficiencyIPNonhfeD[icentr]){
+ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]);
+ weightErr[icentr] = 0;
+ }
+ else{
+ printf("Fit failed on dalitz IP cut efficiency for centrality %d\n",icentr);
+ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1);
+ weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
+ }
+ }
+ }
+ else{
+ if(source==0){
+ if(fEfficiencyIPCharmD[icentr]){
+ weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
+ weightErr[icentr] = 0;
+ }
+ else{
+ printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
+ weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
+ weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
+ }
+ }
+ else if(source==2){
+ if(fBeamType==0){
+ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
+ weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
+ }
+ else{
+ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion
+ weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
+ }
+ }
+ else if(source==3){
+ if(fBeamType==0){
+ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
+ weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
+ }
+ else{
+ weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz
+ weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
+ }
+ }
+ }
+
+ if(fBeamType==1){
+ contents[0]=icentr;
+ contents[1]=pt[0];
+ ibin[0]=icentr;
+ ibin[1]=i+1;
+ }
+ if(fBeamType==0){
+ contents[0]=pt[0];
+ ibin[0]=i+1;
+ }
+
+ pideff->Fill(contents,weight[icentr]);
+ pideff->SetBinError(ibin,weightErr[icentr]);
+ }
+ }
- return ipcut;
+ Int_t nDimSparse = pideff->GetNdimensions();
+ Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
+ Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
+
+ for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+ binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
+ nBins *= binsvar[iVar];
+ }
+
+ Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+ // loop that sets 0 error in each bin
+ for (Long_t iBin=0; iBin<nBins; iBin++) {
+ Long_t bintmp = iBin ;
+ for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+ binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+ bintmp /= binsvar[iVar] ;
+ }
+ }
+
+ delete[] binsvar;
+ delete[] binfill;
+
+
+ return pideff;
}
-//____________________________________________________________________________
-THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
+//__________________________________________________________________________
+AliCFDataGrid *AliHFEspectrum::GetRawBspectra2ndMethod(){
+ //
+ // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
+ //
+ Int_t ptpr = 0;
+ Int_t nDim = 1;
+ if(fBeamType==0)
+ {
+ ptpr=0;
+ }
+ if(fBeamType==1)
+ {
+ ptpr=1;
+ nDim=2;
+ }
+
+ const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
+ const Int_t nCentralitybinning=11;//number of centrality bins
+ Double_t kPtRange[19] = {0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 2, 2.5, 3, 4, 5, 6, 8, 12, 16, 20};
+
+ Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
+ Int_t nBin[1] = {nPtbinning1};
+ Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
+
+ AliCFDataGrid *rawBeautyContainer;
+ if(fBeamType==0) rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
+ else rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBinPbPb);
+ // printf("number of bins= %d\n",bins[0]);
+
+
+
+
+ THnSparseF *brawspectra;
+ if(fBeamType==0) brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
+ else brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBinPbPb);
+ if(fBeamType==0) brawspectra->SetBinEdges(0, kPtRange);
+ if(fBeamType==1)
+ {
+ // brawspectra->SetBinEdges(0, centralityBins);
+ brawspectra->SetBinEdges(0, kCentralityRange);
+ brawspectra->SetBinEdges(1, kPtRange);
+ }
+
+ Double_t pt[1];
+ Double_t yields= 0.;
+ Double_t valuesb[2];
+
+ //Int_t looppt=nBin[0];
+ Int_t loopcentr=1;
+ if(fBeamType==1)
+ {
+ loopcentr=nBinPbPb[0];
+ }
+
+ for(int icentr=0; icentr<loopcentr; icentr++)
+ {
+
+ for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
+ pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
+
+ yields = fBSpectrum2ndMethod->GetBinContent(i+1);
+
+ if(fBeamType==1)
+ {
+ valuesb[0]=icentr;
+ valuesb[1]=pt[0];
+ }
+ if(fBeamType==0) valuesb[0]=pt[0];
+ brawspectra->Fill(valuesb,yields);
+ }
+ }
+
+
+
+ Int_t nDimSparse = brawspectra->GetNdimensions();
+ Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
+ Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
+
+ for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+ binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
+ nBins *= binsvar[iVar];
+ }
+
+ Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+ // loop that sets 0 error in each bin
+ for (Long_t iBin=0; iBin<nBins; iBin++) {
+ Long_t bintmp = iBin ;
+ for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+ binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+ bintmp /= binsvar[iVar] ;
+ }
+ brawspectra->SetBinError(binfill,0.); // put 0 everywhere
+ }
+
+
+ rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
+ TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(ptpr);
+
+ new TCanvas;
+ fBSpectrum2ndMethod->SetMarkerStyle(24);
+ fBSpectrum2ndMethod->Draw("p");
+ hRawBeautySpectra->SetMarkerStyle(25);
+ hRawBeautySpectra->Draw("samep");
+
+ delete[] binfill;
+ delete[] binsvar;
+
+ return rawBeautyContainer;
+}
+
+//__________________________________________________________________________
+void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){
+ //
+ // Calculate non HFE sys
//
- // Return PID x IP cut efficiency
//
- const Int_t nDim=1;
- Int_t nBin[nDim] = {44};
- const Double_t kPtbound[2] = {0.1, 20.};
+ if(!fNonHFEsyst)
+ return;
- Double_t* binEdges[nDim];
- binEdges[0] = AliHFEtools::MakeLogarithmicBinning(nBin[0], kPtbound[0], kPtbound[1]);
+ Double_t evtnorm[1] = {0.0};
+ if(fNMCbgEvents[0]>0) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
+
+ AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels];
+ AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels];
- THnSparseF *pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
- for(Int_t idim = 0; idim < nDim; idim++)
- pideff->SetBinEdges(idim, binEdges[idim]);
+ AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors
+ AliCFDataGrid *bgNonHFEGrid[kBgLevels];
+ AliCFDataGrid *bgConvGrid[kBgLevels];
- Double_t pt[1];
- Double_t weight = 1.0;
- Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
- for(int i=0; i<nBin[0]; i++){
- pt[0]=(binEdges[0][i]+binEdges[0][i+1])/2.;
- Double_t tofeff = 0.9885*(0.8551+0.0070*log(pt[0])-0.0014*log(pt[0])*log(pt[0]))*tanh(0.8884*pt[0]+0.6061); // parameterized TOF PID eff
-
- if(source==0) weight = tofeff*trdtpcPidEfficiency*0.299*(0.307+0.176*log(pt[0])+0.0176*log(pt[0])*log(pt[0]))*tanh(96.3*pt[0]-19.7); // for 3 sigma cut
- //if(source==0) weight = tofeff*trdtpcPidEfficiency*0.477*(0.3757+0.184*log(pt[0])-0.0013*log(pt[0])*log(pt[0]))*tanh(436.013*pt[0]-96.2137); // for 2 sigma cut
- //if(source==2) weight = tofeff*trdtpcPidEfficiency*0.5575*(0.3594-0.3051*log(pt[0])+0.1597*log(pt[0])*log(pt[0]))*tanh(8.2436*pt[0]-0.8125);
- //if(source==3) weight = tofeff*trdtpcPidEfficiency*0.0937*(0.4449-0.3769*log(pt[0])+0.5031*log(pt[0])*log(pt[0]))*tanh(6.4063*pt[0]+3.1425); // multiply ip cut eff for Dalitz/dielectrons
- if(source==2) weight = tofeff*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1); // multiply ip cut eff for conversion electrons
- if(source==3) weight = tofeff*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // multiply ip cut eff for Dalitz/dielectrons
- printf("tofeff= %lf trdtpcPidEfficiency= %lf conversion= %lf nonhfe= %lf\n",tofeff,trdtpcPidEfficiency,fConversionEff->GetBinContent(i+1),fNonHFEEff->GetBinContent(i+1));
-
- pideff->Fill(pt,weight);
- }
- Int_t* ibins[nDim];
- for(Int_t ivar = 0; ivar < nDim; ivar++)
- ibins[ivar] = new Int_t[nBin[ivar] + 1];
- pideff->SetBinError(ibins[0],0);
+ Int_t stepbackground = 3;
+ Int_t* bins=new Int_t[1];
+ const Char_t *bgBase[2] = {"pi0","eta"};
+
+ bins[0]=fConversionEff[centrality]->GetNbinsX();
+
+ AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
+ AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
- return pideff;
+ for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
+ for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
+ convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),*fConvSourceContainer[iSource][iLevel][centrality],stepbackground);
+ weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
+ convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
+
+ nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),*fNonHFESourceContainer[iSource][iLevel][centrality],stepbackground);
+ weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
+ nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
+ }
+
+ bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
+ for(Int_t iSource = 2; iSource < kElecBgSources; iSource++){
+ bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
+ }
+ if(!fEtaSyst)
+ bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]);
+
+ bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
+ for(Int_t iSource = 2; iSource < kElecBgSources; iSource++){//add other sources to pi0, to get overall background from all meson decays, exception: eta (independent error calculation)
+ bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
+ }
+ if(!fEtaSyst)
+ bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]);
+
+ bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
+ bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]);
+ if(fEtaSyst){
+ bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source
+ bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]);
+ }
+ }
+
+
+ //Now subtract the mean from upper, and lower from mean container to get the error based on the pion yield uncertainty (-> this error sums linearly, since its contribution to all meson yields is correlated; exception: eta errors in pp 7 TeV sum with others the gaussian way, as they are independent from pi0)
+ AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper
+ TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper
+ for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base
+ bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone();
+ bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.);
+ bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone();
+ bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.);
+
+ //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate
+
+ hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0);
+ hBaseErrors[iErr][0]->Scale(-1.);
+ hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr]));
+ hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0);
+ if(!fEtaSyst)break;
+ }
+
+
+ //Calculate the scaling errors for electrons from all mesons except for pions (and in pp 7 TeV case eta): square sum of (0.3 * best yield estimate), where 0.3 is the error generally assumed for m_t scaling
+ TH1D *hSpeciesErrors[kElecBgSources-1];
+ for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
+ if(fEtaSyst && (iSource == 1))continue;
+ hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
+ TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
+ hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
+ hSpeciesErrors[iSource-1]->Scale(0.3);
+ }
+
+ //Int_t firstBgSource = 0;//if eta systematics are not from scaling
+ //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined!
+ TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone();
+ TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone();
+ TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone();
+
+ TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
+ TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
+
+ for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
+ Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp;
+ pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin);
+ pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin);
+ if(fEtaSyst){
+ etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin);
+ etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin);
+ }
+ else{ etaErrLow = etaErrUp = 0.;}
+
+ Double_t sqrsumErrs= 0;
+ for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
+ if(fEtaSyst && (iSource == 1))continue;
+ Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
+ sqrsumErrs+=(scalingErr*scalingErr);
+ }
+ for(Int_t iErr = 0; iErr < 2; iErr++){
+ for(Int_t iLevel = 0; iLevel < 2; iLevel++){
+ hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin));
+ }
+ if(!fEtaSyst)break;
+ }
+ hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs));
+ hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs));
+ hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
+
+ hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
+ hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));
+ }
+
+
+ TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
+ cPiErrors->cd();
+ cPiErrors->SetLogx();
+ cPiErrors->SetLogy();
+ hBaseErrors[0][0]->Draw();
+ //hBaseErrors[0][1]->SetMarkerColor(kBlack);
+ //hBaseErrors[0][1]->SetLineColor(kBlack);
+ //hBaseErrors[0][1]->Draw("SAME");
+ if(fEtaSyst){
+ hBaseErrors[1][0]->Draw("SAME");
+ hBaseErrors[1][0]->SetMarkerColor(kBlack);
+ hBaseErrors[1][0]->SetLineColor(kBlack);
+ //hBaseErrors[1][1]->SetMarkerColor(13);
+ //hBaseErrors[1][1]->SetLineColor(13);
+ //hBaseErrors[1][1]->Draw("SAME");
+ }
+ //hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
+ //hOverallBinScaledErrsUp->SetLineColor(kBlue);
+ //hOverallBinScaledErrsUp->Draw("SAME");
+ hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
+ hOverallBinScaledErrsLow->SetLineColor(kGreen);
+ hOverallBinScaledErrsLow->Draw("SAME");
+ hScalingErrors->SetLineColor(kBlue);
+ hScalingErrors->Draw("SAME");
+
+ TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
+ lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error");
+ //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error");
+ if(fEtaSyst){
+ lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error");
+ //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error");
+ }
+ lPiErr->AddEntry(hScalingErrors, "scaling error");
+ lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
+ //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
+ lPiErr->Draw("SAME");
+
+ //Normalize errors
+ TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
+ TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
+ hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
+ hLowSystScaled->Scale(evtnorm[0]);
+ TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
+ TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
+ //histograms to be normalized to TGraphErrors
+ CorrectFromTheWidth(hNormAllSystErrUp);
+ CorrectFromTheWidth(hNormAllSystErrLow);
+
+ TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
+ cNormOvErrs->cd();
+ cNormOvErrs->SetLogx();
+ cNormOvErrs->SetLogy();
+
+ TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
+ TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
+ gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
+ gOverallSystErrUp->SetMarkerColor(kBlack);
+ gOverallSystErrUp->SetLineColor(kBlack);
+ gOverallSystErrLow->SetMarkerColor(kRed);
+ gOverallSystErrLow->SetLineColor(kRed);
+ gOverallSystErrUp->Draw("AP");
+ gOverallSystErrLow->Draw("PSAME");
+ TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
+ lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
+ lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
+ lAllSys->Draw("same");
+
+
+ AliCFDataGrid *bgYieldGrid;
+ if(fEtaSyst){
+ bgLevelGrid[0][0]->Add(bgLevelGrid[1][0]);//Addition of the eta background best estimate to the rest. Needed to be separated for error treatment - now overall background necessary! If no separate eta systematics exist, the corresponding grid has already been added before.
+ }
+ bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone();
+
+ TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
+ TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
+ hRelErrUp->Divide(hBgYield);
+ TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
+ hRelErrLow->Divide(hBgYield);
+
+ TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
+ cRelErrs->cd();
+ cRelErrs->SetLogx();
+ hRelErrUp->SetTitle("Relative error of non-HFE background yield");
+ hRelErrUp->Draw();
+ hRelErrLow->SetLineColor(kBlack);
+ hRelErrLow->Draw("SAME");
+
+ TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
+ lRel->AddEntry(hRelErrUp, "upper");
+ lRel->AddEntry(hRelErrLow, "lower");
+ lRel->Draw("SAME");
+
+ //CorrectFromTheWidth(hBgYield);
+ //hBgYield->Scale(evtnorm[0]);
+
+
+ //write histograms/TGraphs to file
+ TFile *output = new TFile("systHists.root","recreate");
+
+ hBgYield->SetName("hBgYield");
+ hBgYield->Write();
+ hRelErrUp->SetName("hRelErrUp");
+ hRelErrUp->Write();
+ hRelErrLow->SetName("hRelErrLow");
+ hRelErrLow->Write();
+ hUpSystScaled->SetName("hOverallSystErrUp");
+ hUpSystScaled->Write();
+ hLowSystScaled->SetName("hOverallSystErrLow");
+ hLowSystScaled->Write();
+ gOverallSystErrUp->SetName("gOverallSystErrUp");
+ gOverallSystErrUp->Write();
+ gOverallSystErrLow->SetName("gOverallSystErrLow");
+ gOverallSystErrLow->Write();
+
+ output->Close();
+ delete output;
+
+}
+
+//____________________________________________________________
+void AliHFEspectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){
+
+ //
+ // Unfold backgrounds to check its sanity
+ //
+
+ AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
+ //AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
+ if(!mcContainer){
+ AliError("MC Container not available");
+ return;
+ }
+
+ if(!fCorrelation){
+ AliError("No Correlation map available");
+ return;
+ }
+
+ // Data
+ AliCFDataGrid *dataGrid = 0x0;
+ dataGrid = bgsubpectrum;
+
+ // Guessed
+ AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
+ THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
+
+ // Efficiency
+ AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
+ efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue);
+
+ // Unfold background spectra
+ Int_t nDim=1;
+ if(fBeamType==0)nDim = 1;
+ if(fBeamType==1)nDim = 2;
+ AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
+ if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
+ unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
+ if(fSetSmoothing) unfolding.UseSmoothing();
+ unfolding.Unfold();
+
+ // Results
+ THnSparse* result = unfolding.GetUnfolded();
+ TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600);
+ if(fBeamType==1)
+ {
+ ctest->Divide(2,2);
+ ctest->cd(1);
+ result->GetAxis(0)->SetRange(1,1);
+ TH1D* htest1=(TH1D*)result->Projection(0);
+ htest1->Draw();
+ ctest->cd(2);
+ result->GetAxis(0)->SetRange(1,1);
+ TH1D* htest2=(TH1D*)result->Projection(1);
+ htest2->Draw();
+ ctest->cd(3);
+ result->GetAxis(0)->SetRange(6,6);
+ TH1D* htest3=(TH1D*)result->Projection(0);
+ htest3->Draw();
+ ctest->cd(4);
+ result->GetAxis(0)->SetRange(6,6);
+ TH1D* htest4=(TH1D*)result->Projection(1);
+ htest4->Draw();
+
+ }
+
+
+
+
+
+ TGraphErrors* unfoldedbgspectrumD = Normalize(result);
+ if(!unfoldedbgspectrumD) {
+ AliError("Unfolded background spectrum doesn't exist");
+ }
+ else{
+ TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate");
+ if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum");
+
+ if(fBeamType==1)
+ {
+ Int_t centr=1;
+ result->GetAxis(0)->SetRange(centr,centr);
+ unfoldedbgspectrumD = Normalize(result,centr-1);
+ unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr0_20");
+ centr=6;
+ result->GetAxis(0)->SetRange(centr,centr);
+ unfoldedbgspectrumD = Normalize(result,centr-1);
+ unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr40_80");
+ }
+
+ file->Close();
+ }
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff