/************************************************************************** * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /* $Id: $ */ ///////////////////////////////////////////////////////////// // class to average D meson -hadron correlations // // Author: A. Rossi, andrea.rossi@cern.ch ///////////////////////////////////////////////////////////// #include "AliHFDmesonCorrAverage.h" #include "AliHFDhadronCorrSystUnc.h" #include #include #include #include //using std::cout; //using std::endl; ClassImp(AliHFDmesonCorrAverage) AliHFDmesonCorrAverage::AliHFDmesonCorrAverage() : TNamed(), fSystDzero(0), fSystDstar(0), fSystDplus(0), fSystDaverage(0), fincludeDzero(kFALSE), fincludeDstar(kFALSE), fincludeDplus(kFALSE), fmethod(), fptminD(), fptmaxD(), fptminAsso(), fptmaxAsso(), fhDzero(0x0), fhDstar(0x0), fhDplus(0x0), fhDaverage(0x0), fgrTotSystAverage(0x0), fgrFDSystAverage(0x0), fgrNonFDSystAverage(0x0), fweightsDzeroStat(0x0), fweightsDstarStat(0x0), fweightsDplusStat(0x0), fweightsDzeroSystYield(0x0), fweightsDstarSystYield(0x0), fweightsDplusSystYield(0x0), fweightsDzeroSystBkg(0x0), fweightsDstarSystBkg(0x0), fweightsDplusSystBkg(0x0), fnbinsphi(), fsys(-2), fyear(-2), fSystAlreadySet(kFALSE) {// standard constructor } AliHFDmesonCorrAverage::AliHFDmesonCorrAverage(const char* name) : TNamed(name,name), fSystDzero(0), fSystDstar(0), fSystDplus(0), fSystDaverage(0), fincludeDzero(kFALSE), fincludeDstar(kFALSE), fincludeDplus(kFALSE), fmethod(), fptminD(), fptmaxD(), fptminAsso(), fptmaxAsso(), fhDzero(0x0), fhDstar(0x0), fhDplus(0x0), fhDaverage(0x0), fgrTotSystAverage(0x0), fgrFDSystAverage(0x0), fgrNonFDSystAverage(0x0), fweightsDzeroStat(0x0), fweightsDstarStat(0x0), fweightsDplusStat(0x0), fweightsDzeroSystYield(0x0), fweightsDstarSystYield(0x0), fweightsDplusSystYield(0x0), fweightsDzeroSystBkg(0x0), fweightsDstarSystBkg(0x0), fweightsDplusSystBkg(0x0), fnbinsphi(), fsys(-2), fyear(-2), fSystAlreadySet(kFALSE) {// default constructor } AliHFDmesonCorrAverage::~AliHFDmesonCorrAverage(){ delete fSystDzero; delete fSystDstar; delete fSystDplus; delete fSystDaverage; delete fhDzero; delete fhDstar; delete fhDplus; delete fhDaverage; delete fgrTotSystAverage; delete fgrFDSystAverage; delete fgrNonFDSystAverage; delete fweightsDzeroStat; delete fweightsDstarStat; delete fweightsDplusStat; delete fweightsDzeroSystYield; delete fweightsDstarSystYield; delete fweightsDplusSystYield; delete fweightsDzeroSystBkg; delete fweightsDstarSystBkg; delete fweightsDplusSystBkg; } Bool_t AliHFDmesonCorrAverage::InitSystematicUncertainty(Int_t system,Int_t year){ if(!fSystAlreadySet){ if(system!=-1){ if(system==0){ //pp if(fincludeDzero){ if(year==2010){ fSystDzero=new AliHFDhadronCorrSystUnc("fSystDzero"); fSystDzero->InitStandardUncertaintiesPP2010(0,(fptminD+fptmaxD)*0.5,fptminAsso); fSystDzero->BuildTotalUncHisto(); fSystDzero->BuildTotalNonFDUncHisto(); fSystDzero->BuildGraphsUnc(fhDzero); fSystDzero->BuildGraphsRelUnc(); } else { Printf("No values for syst unc foreseen for this dataset"); } } if(fincludeDstar){ if(year==2010){ fSystDstar=new AliHFDhadronCorrSystUnc("fSystDstar"); fSystDstar->InitStandardUncertaintiesPP2010(1,(fptminD+fptmaxD)*0.5,fptminAsso); fSystDstar->BuildTotalUncHisto(); fSystDstar->BuildTotalNonFDUncHisto(); fSystDstar->BuildGraphsUnc(fhDstar); fSystDstar->BuildGraphsRelUnc(); } else { Printf("No values for syst unc foreseen for this dataset"); } } if(fincludeDplus){ if(year==2010){ fSystDplus=new AliHFDhadronCorrSystUnc("fSystDplus"); fSystDplus->InitStandardUncertaintiesPP2010(2,(fptminD+fptmaxD)*0.5,fptminAsso); fSystDplus->BuildTotalUncHisto(); fSystDplus->BuildTotalNonFDUncHisto(); fSystDplus->BuildGraphsUnc(fhDplus); fSystDplus->BuildGraphsRelUnc(); } else { Printf("No values for syst unc foreseen for this dataset"); } } } else if(system==1){ //p-Pb if(fincludeDzero){ if(year==2013){ fSystDzero=new AliHFDhadronCorrSystUnc("fSystDzero"); fSystDzero->InitStandardUncertaintiesPPb2013(0,(fptminD+fptmaxD)*0.5,fptminAsso); fSystDzero->BuildTotalUncHisto(); fSystDzero->BuildTotalNonFDUncHisto(); fSystDzero->BuildGraphsUnc(fhDzero); fSystDzero->BuildGraphsRelUnc(); } else { Printf("No values for syst unc foreseen for this dataset"); } } if(fincludeDstar){ if(year==2013){ fSystDstar=new AliHFDhadronCorrSystUnc("fSystDstar"); fSystDstar->InitStandardUncertaintiesPPb2013(1,(fptminD+fptmaxD)*0.5,fptminAsso); fSystDstar->BuildTotalUncHisto(); fSystDstar->BuildTotalNonFDUncHisto(); fSystDstar->BuildGraphsUnc(fhDstar); fSystDstar->BuildGraphsRelUnc(); } else { Printf("No values for syst unc foreseen for this dataset"); } } if(fincludeDplus){ if(year==2013){ fSystDplus=new AliHFDhadronCorrSystUnc("fSystDplus"); fSystDplus->InitStandardUncertaintiesPPb2013(2,(fptminD+fptmaxD)*0.5,fptminAsso); fSystDplus->BuildTotalUncHisto(); fSystDplus->BuildTotalNonFDUncHisto(); fSystDplus->BuildGraphsUnc(fhDplus); fSystDplus->BuildGraphsRelUnc(); } else { Printf("No values for syst unc foreseen for this dataset"); } } } else { Printf("Cannot initiate syst uncertainties: wrong system selected"); return kFALSE; } } } fSystDaverage=new AliHFDhadronCorrSystUnc("fSystDaverage"); if(fincludeDzero)fSystDaverage->SetHistoTemplate(fSystDzero->GetHistoTemplate(),"fhDeltaPhiTemplateDaverage"); else if(fincludeDstar)fSystDaverage->SetHistoTemplate(fSystDstar->GetHistoTemplate(),"fhDeltaPhiTemplateDaverage"); else if(fincludeDplus)fSystDaverage->SetHistoTemplate(fSystDplus->GetHistoTemplate(),"fhDeltaPhiTemplateDaverage"); fSystDaverage->InitEmptyHistosFromTemplate(); return kTRUE; } void AliHFDmesonCorrAverage::InitAverageHisto(TH1D *h){ fhDaverage=(TH1D*)h->Clone("fhDaverage"); fhDaverage->SetXTitle("#Delta#varphi = #varphi_{assoc} - #varphi_{D}"); fhDaverage->SetYTitle("#frac{1}{N_{D}}#frac{dN^{h}}{d#Delta#varphi} (rad^{-1})"); fhDaverage->SetTitle(""); fhDaverage->SetMinimum(0); fhDaverage->Reset(0); fhDaverage->SetLineWidth(2); fhDaverage->SetLineColor(kBlack); fhDaverage->SetMarkerStyle(20); fhDaverage->SetMarkerSize(1.2); fhDaverage->SetMarkerColor(kBlack); } void AliHFDmesonCorrAverage::CalculateAverage(){ if(fmethod!=10&&fmethod!=20){ Printf("This option for the calculation of the average has not been implemented yet"); return; } // check that histos were set if(fincludeDzero&&(!fhDzero)){ Printf("Dzero histo not set"); return; } if(fincludeDstar&&(!fhDstar)){ Printf("Dstar histo not set"); return; } if(fincludeDplus&&(!fhDplus)){ Printf("Dplus histo not set"); return; } // check that they have the same binning if(fincludeDzero&&fincludeDstar){ if(fhDzero->GetNbinsX()!=fhDstar->GetNbinsX()){ Printf("Dzero and Dstar histos have different binning"); return ; } fnbinsphi=fhDzero->GetNbinsX(); for(Int_t j=1;j<=fnbinsphi;j++){ if(TMath::Abs(fhDzero->GetBinLowEdge(j)-fhDstar->GetBinLowEdge(j))>0.001){ Printf("Dzero and Dstar histos have different binning"); return; } } InitAverageHisto(fhDzero); } if(fincludeDzero&&fincludeDplus){ if(fhDzero->GetNbinsX()!=fhDplus->GetNbinsX()){ Printf("Dzero and Dplus histos have different binning"); return ; } fnbinsphi=fhDzero->GetNbinsX(); for(Int_t j=1;j<=fnbinsphi;j++){ if(TMath::Abs(fhDzero->GetBinLowEdge(j)-fhDplus->GetBinLowEdge(j))>0.001){ Printf("Dzero and Dplus histos have different binning"); return; } } if(!fhDaverage)InitAverageHisto(fhDzero); } if(fincludeDstar&&fincludeDplus){ if(fhDstar->GetNbinsX()!=fhDplus->GetNbinsX()){ Printf("Dstar and Dplus histos have different binning"); return ; } fnbinsphi=fhDstar->GetNbinsX(); for(Int_t j=1;j<=fnbinsphi;j++){ if(TMath::Abs(fhDstar->GetBinLowEdge(j)-fhDplus->GetBinLowEdge(j))>0.001){ Printf("Dplus and Dstar histos have different binning"); return; } } if(!fhDaverage)InitAverageHisto(fhDstar); } if(!(fincludeDstar||fincludeDplus||fincludeDzero)){ Printf("No mesons choosen for average"); return ; } Int_t nmeson=0; if(fincludeDzero)nmeson++; if(fincludeDstar)nmeson++; if(fincludeDplus)nmeson++; printf("Settin syst \n"); if(!InitSystematicUncertainty(fsys,fyear)){ Printf("Cannot init syst uncertainties \n"); return; } SetWeights(); printf("Weights set \n"); for(Int_t j=1;j<=fnbinsphi;j++){ Double_t value=0.,errStatValue=0.,errSystYieldValue=0.,errSystBkgValue=0.,weight=0.; Double_t systMCclosureMin=0.,systMCcorrectionsMin=0.,systFDmin=0.,systMCDefficiencyMin=0.,systSecContaminationMin=0.; Double_t systMCclosureMax=0.,systMCcorrectionsMax=0.,systFDmax=0.,systMCDefficiencyMax=0.,systSecContaminationMax=0.; Double_t sumweights=0; Double_t errsyst; if(fincludeDzero){ // stat error + yield unc + bkg subtr weight=1./(1./fweightsDzeroStat[j-1]+1./fweightsDzeroSystYield[j-1]+1./fweightsDzeroSystBkg[j-1]);// need to do this way since we stored separately the stat and syst weight (=1/variance) value+=fhDzero->GetBinContent(j)*weight; errStatValue+=1./fweightsDzeroStat[j-1]*weight*weight; errSystYieldValue+=1./fweightsDzeroSystYield[j-1]*weight*weight; errSystBkgValue+=1./fweightsDzeroSystBkg[j-1]*weight*weight; sumweights+=weight; printf("Dzero the value is: %f, weight: %f \n",value, weight); // MCcorrections : associated tracks, correlated errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCcorrectionsMin()->GetBinContent(j)); systMCcorrectionsMin+=errsyst*weight; errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCcorrectionsMax()->GetBinContent(j); systMCcorrectionsMax+=errsyst*weight; printf(" Dzero trackeff the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCcorrectionsMin()->GetBinContent(j), systMCcorrectionsMin); printf(" Dzero trackeff the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCcorrectionsMax()->GetBinContent(j), systMCcorrectionsMax); // MCDefficiency : uncorrelated errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCDefficiencyMin()->GetBinContent(j); systMCDefficiencyMin+=errsyst*errsyst*weight*weight; errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCDefficiencyMax()->GetBinContent(j); systMCDefficiencyMax+=errsyst*errsyst*weight*weight; // SecContamination : correlated errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoSecContaminationMin()->GetBinContent(j)); systSecContaminationMin+=errsyst*weight; errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoSecContaminationMax()->GetBinContent(j); systSecContaminationMax+=errsyst*weight; // MC closure: fully correlated errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCclosureTestMin()->GetBinContent(j)); systMCclosureMin+=errsyst*weight; errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCclosureTestMax()->GetBinContent(j); systMCclosureMax+=errsyst*weight; printf(" Dzero Mcclosure the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCclosureTestMin()->GetBinContent(j), systMCclosureMin); printf(" Dzero Mcclosure the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCclosureTestMax()->GetBinContent(j), systMCclosureMax); // FD (assume full correlations) errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoFDmin()->GetBinContent(j)); systFDmin+=errsyst*weight; errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoFDmax()->GetBinContent(j); systFDmax+=errsyst*weight; printf(" Dzero feeddown the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoFDmin()->GetBinContent(j), systFDmin); printf(" Dzero feeddown the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoFDmax()->GetBinContent(j), systFDmax); printf("Dzero the value is: %f, weight: %f \n",value, weight); } if(fincludeDstar){ // stat error + yield unc + bkg subtr weight=1./(1./fweightsDstarStat[j-1]+1./fweightsDstarSystYield[j-1]+1./fweightsDstarSystBkg[j-1]);// need to do this way since we stored separately the stat and syst weight (=1/variance) value+=fhDstar->GetBinContent(j)*weight; errStatValue+=1./fweightsDstarStat[j-1]*weight*weight; errSystYieldValue+=1./fweightsDstarSystYield[j-1]*weight*weight; errSystBkgValue+=1./fweightsDstarSystBkg[j-1]*weight*weight; sumweights+=weight; printf("Dstar the value is: %f, weight: %f \n",value, weight); // MCcorrections: associated tracks, correlated errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCcorrectionsMin()->GetBinContent(j)); systMCcorrectionsMin+=errsyst*weight; errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCcorrectionsMax()->GetBinContent(j)); systMCcorrectionsMax+=errsyst*weight; printf(" Dstar trackeff the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCcorrectionsMin()->GetBinContent(j), systMCcorrectionsMin); printf(" Dstar trackeff the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCcorrectionsMax()->GetBinContent(j), systMCcorrectionsMax); // MCDefficiency : uncorrelated errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCDefficiencyMin()->GetBinContent(j); systMCDefficiencyMin+=errsyst*errsyst*weight*weight; errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCDefficiencyMax()->GetBinContent(j); systMCDefficiencyMax+=errsyst*errsyst*weight*weight; // SecContamination : correlated errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoSecContaminationMin()->GetBinContent(j)); systSecContaminationMin+=errsyst*weight; errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoSecContaminationMax()->GetBinContent(j); systSecContaminationMax+=errsyst*weight; // MC closure errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCclosureTestMin()->GetBinContent(j)); systMCclosureMin+=errsyst*weight; errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCclosureTestMax()->GetBinContent(j); systMCclosureMax+=errsyst*weight; printf(" Dstar Mcclosure the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCclosureTestMin()->GetBinContent(j), systMCclosureMin); printf(" Dstar Mcclosure the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCclosureTestMax()->GetBinContent(j), systMCclosureMax); // FD (assume full correlations) errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoFDmin()->GetBinContent(j)); systFDmin+=errsyst*weight; errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoFDmax()->GetBinContent(j); systFDmax+=errsyst*weight; printf(" Dstar feeddown the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoFDmin()->GetBinContent(j), systFDmin); printf(" Dstar feeddown the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoFDmax()->GetBinContent(j), systFDmax); printf("Dstar the value is: %f, weight: %f \n",value, weight); } if(fincludeDplus){ // stat error + yield unc + bkg subtr weight=1./(1./fweightsDplusStat[j-1]+1./fweightsDplusSystYield[j-1]+1./fweightsDplusSystBkg[j-1]);// need to do this way since we stored separately the stat and syst weight (=1/variance) value+=fhDplus->GetBinContent(j)*weight; errStatValue+=1./fweightsDplusStat[j-1]*weight*weight; errSystYieldValue+=1./fweightsDplusSystYield[j-1]*weight*weight; errSystBkgValue+=1./fweightsDplusSystBkg[j-1]*weight*weight; sumweights+=weight; printf("Dplus the value is: %f, weight: %f \n",value, weight); // MCcorrections : associated tracsk, correlated errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCcorrectionsMin()->GetBinContent(j)); systMCcorrectionsMin+=errsyst*weight; errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCcorrectionsMax()->GetBinContent(j); systMCcorrectionsMax+=errsyst*weight; printf(" Dplus trackeff the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCcorrectionsMin()->GetBinContent(j), systMCcorrectionsMin); printf(" Dplus trackeff the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCcorrectionsMax()->GetBinContent(j), systMCcorrectionsMax); // MCDefficiency : uncorrelated errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCDefficiencyMin()->GetBinContent(j); systMCDefficiencyMin+=errsyst*errsyst*weight*weight; errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCDefficiencyMax()->GetBinContent(j); systMCDefficiencyMax+=errsyst*errsyst*weight*weight; printf(" Dplus cutvariat the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCDefficiencyMin()->GetBinContent(j), systMCDefficiencyMin); printf(" Dplus cutvariat the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCDefficiencyMax()->GetBinContent(j), systMCDefficiencyMax); // SecContamination : correlated errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoSecContaminationMin()->GetBinContent(j)); systSecContaminationMin+=errsyst*weight; errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoSecContaminationMax()->GetBinContent(j); systSecContaminationMax+=errsyst*weight; // MC closure errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCclosureTestMin()->GetBinContent(j)); systMCclosureMin+=errsyst*weight; errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCclosureTestMax()->GetBinContent(j); systMCclosureMax+=errsyst*weight; printf(" Dplus Mcclosure the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCclosureTestMin()->GetBinContent(j), systMCclosureMin); printf(" Dplus Mcclosure the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCclosureTestMax()->GetBinContent(j), systMCclosureMax); // FD (assume full correlations) errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoFDmin()->GetBinContent(j)); systFDmin+=errsyst*weight; errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoFDmax()->GetBinContent(j); systFDmax+=errsyst*weight; printf(" Dplus feeddown the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoFDmin()->GetBinContent(j), systMCclosureMin); printf(" Dplus feeddown the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoFDmax()->GetBinContent(j), systMCclosureMax); printf("Dplus the value is: %f, weight: %f \n",value, weight); } value/=sumweights; errStatValue/=sumweights*sumweights;// was: ((Double_t)nmeson*sumweights); errSystYieldValue/=sumweights*sumweights;// was: ((Double_t)nmeson*sumweights); errSystBkgValue/=sumweights*sumweights; // was: ((Double_t)nmeson*sumweights); fhDaverage->SetBinContent(j,value); fhDaverage->SetBinError(j,TMath::Sqrt(errStatValue)); // Settting unc on yield and back sub fSystDaverage->GetHistoYieldUnc()->SetBinContent(j,TMath::Sqrt(errSystYieldValue)/TMath::Abs(value)); fSystDaverage->GetHistoBackSubUncMin()->SetBinContent(j,-TMath::Sqrt(errSystBkgValue)/TMath::Abs(value)); fSystDaverage->GetHistoBackSubUncMax()->SetBinContent(j,TMath::Sqrt(errSystBkgValue)/TMath::Abs(value)); // MCcorrections, associated tracks if(systMCcorrectionsMin>=0.){ systMCcorrectionsMin=systMCcorrectionsMin/sumweights; } else { Printf("Negative sum in error calculation \n"); return; } if(systMCcorrectionsMax>=0.){ systMCcorrectionsMax=systMCcorrectionsMax/sumweights; } else { Printf("Negative sum in error calculation \n"); return; } // MCDefficiency if(systMCDefficiencyMin>=0.){ systMCDefficiencyMin=TMath::Sqrt(systMCDefficiencyMin)/sumweights; } else { Printf("Negative sum in error calculation \n"); return; } if(systMCDefficiencyMax>=0.){ systMCDefficiencyMax=TMath::Sqrt(systMCDefficiencyMax)/sumweights; } else { Printf("Negative sum in error calculation \n"); return; } // SecContamination if(systSecContaminationMin>=0.){ systSecContaminationMin=systSecContaminationMin/sumweights; } else { Printf("Negative sum in error calculation \n"); return; } if(systSecContaminationMax>=0.){ systSecContaminationMax=systSecContaminationMax/sumweights; } else { Printf("Negative sum in error calculation \n"); return; } // MCclosureTest if(systMCclosureMin>=0.){ systMCclosureMin=systMCclosureMin/sumweights; } else { Printf("Negative sum in error calculation \n"); return; } if(systMCclosureMax>=0.){ systMCclosureMax=systMCclosureMax/sumweights; } else { Printf("Negative sum in error calculation \n"); return; } // Feed down if(systFDmin>=0.){ systFDmin=systFDmin/sumweights; } else { Printf("Negative sum in error calculation \n"); return; } if(systFDmax>=0.){ systFDmax=systFDmax/sumweights; } else { Printf("Negative sum in error calculation \n"); return; } value=TMath::Abs(value);// protection to avoid flipping of min syst histo sign fSystDaverage->GetHistoMCcorrectionsMin()->SetBinContent(j,-systMCcorrectionsMin/value); fSystDaverage->GetHistoMCcorrectionsMax()->SetBinContent(j,systMCcorrectionsMax/value); fSystDaverage->GetHistoMCDefficiencyMin()->SetBinContent(j,-systMCDefficiencyMin/value); fSystDaverage->GetHistoMCDefficiencyMax()->SetBinContent(j,systMCDefficiencyMax/value); fSystDaverage->GetHistoSecContaminationMin()->SetBinContent(j,-systSecContaminationMin/value); fSystDaverage->GetHistoSecContaminationMax()->SetBinContent(j,systSecContaminationMax/value); fSystDaverage->GetHistoMCclosureTestMin()->SetBinContent(j,-systMCclosureMin/value); fSystDaverage->GetHistoMCclosureTestMax()->SetBinContent(j,systMCclosureMax/value); fSystDaverage->GetHistoFDmin()->SetBinContent(j,-systFDmin/value); fSystDaverage->GetHistoFDmax()->SetBinContent(j,systFDmax/value); } fSystDaverage->BuildTotalUncHisto(); fSystDaverage->BuildTotalNonFDUncHisto(); fSystDaverage->BuildGraphsUnc(fhDaverage); fSystDaverage->BuildGraphsRelUnc(); return ; } TH1D* AliHFDmesonCorrAverage::ReflectHisto(TH1D *h){ TH1D *h2=new TH1D(Form("%sReflected",h->GetName()),Form("%sReflected",h->GetName()),h->GetNbinsX()/2.,0.,TMath::Pi()); for(Int_t j=1;j<=h->GetNbinsX();j++){ Double_t x=h->GetBinCenter(j); Double_t y0=h->GetBinContent(j); Double_t ey0=h->GetBinError(j); Int_t j2; if(x>0&&xFindBin(x); else if(x<0)j2=h2->FindBin(-1.*x); else if(x>TMath::Pi())j2=h2->FindBin(2.*TMath::Pi()-x); else printf("Point %d excluded \n",j); Double_t y=h2->GetBinContent(j2); Double_t ey=h2->GetBinError(j2); h2->SetBinContent(j2,y+y0); h2->SetBinError(j2,TMath::Sqrt(ey0*ey0+ey*ey)); } return h2; } void AliHFDmesonCorrAverage::SetWeights(){ if(fincludeDzero){ TH1D *hYieldUnc=fSystDzero->GetHistoYieldUnc(); TH1D *hBkgUncMin=fSystDzero->GetHistoBackSubUncMin(); TH1D *hBkgUncMax=fSystDzero->GetHistoBackSubUncMax(); if(fweightsDzeroStat)delete fweightsDzeroStat; if(fweightsDzeroSystYield)delete fweightsDzeroSystYield; if(fweightsDzeroSystBkg)delete fweightsDzeroSystBkg; fweightsDzeroStat=new Double_t[fnbinsphi]; fweightsDzeroSystYield=new Double_t[fnbinsphi]; fweightsDzeroSystBkg=new Double_t[fnbinsphi]; for(Int_t j=0;jGetBinError(j+1)*fhDzero->GetBinError(j+1)); fweightsDzeroSystYield[j]=1./(hYieldUnc->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)*hYieldUnc->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)); //for asymmetric uncertainties... if(TMath::Abs(hBkgUncMin->GetBinContent(j+1)) > TMath::Abs(hBkgUncMax->GetBinContent(j+1))) fweightsDzeroSystBkg[j]=1./(hBkgUncMin->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)*hBkgUncMin->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)); else fweightsDzeroSystBkg[j]=1./(hBkgUncMax->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)*hBkgUncMax->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)); } else{ Printf("This option for the calculation of the average has not been implemented yet"); return; } } } if(fincludeDstar){ TH1D *hYieldUnc=fSystDstar->GetHistoYieldUnc(); TH1D *hBkgUncMin=fSystDstar->GetHistoBackSubUncMin(); TH1D *hBkgUncMax=fSystDstar->GetHistoBackSubUncMax(); if(fweightsDstarStat)delete fweightsDstarStat; if(fweightsDstarSystYield)delete fweightsDstarSystYield; if(fweightsDstarSystBkg)delete fweightsDstarSystBkg; fweightsDstarStat=new Double_t[fnbinsphi]; fweightsDstarSystYield=new Double_t[fnbinsphi]; fweightsDstarSystBkg=new Double_t[fnbinsphi]; for(Int_t j=0;jGetBinError(j+1)*fhDstar->GetBinError(j+1)); fweightsDstarSystYield[j]=1./(hYieldUnc->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)*hYieldUnc->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)); //for asymmetric uncertainties... if(TMath::Abs(hBkgUncMin->GetBinContent(j+1)) > TMath::Abs(hBkgUncMax->GetBinContent(j+1))) fweightsDstarSystBkg[j]=1./(hBkgUncMin->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)*hBkgUncMin->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)); else fweightsDstarSystBkg[j]=1./(hBkgUncMax->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)*hBkgUncMax->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)); } else{ Printf("This option for the calculation of the average has not been implemented yet"); return; } } } if(fincludeDplus){ TH1D *hYieldUnc=fSystDplus->GetHistoYieldUnc(); TH1D *hBkgUncMin=fSystDplus->GetHistoBackSubUncMin(); TH1D *hBkgUncMax=fSystDplus->GetHistoBackSubUncMax(); if(fweightsDplusStat)delete fweightsDplusStat; if(fweightsDplusSystYield)delete fweightsDplusSystYield; if(fweightsDplusSystBkg)delete fweightsDplusSystBkg; fweightsDplusStat=new Double_t[fnbinsphi]; fweightsDplusSystYield=new Double_t[fnbinsphi]; fweightsDplusSystBkg=new Double_t[fnbinsphi]; for(Int_t j=0;jGetBinError(j+1)*fhDplus->GetBinError(j+1)); fweightsDplusSystYield[j]=1./(hYieldUnc->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)*hYieldUnc->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)); //for asymmetric uncertainties... if(TMath::Abs(hBkgUncMin->GetBinContent(j+1)) > TMath::Abs(hBkgUncMax->GetBinContent(j+1))) fweightsDplusSystBkg[j]=1./(hBkgUncMin->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)*hBkgUncMin->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)); else fweightsDplusSystBkg[j]=1./(hBkgUncMax->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)*hBkgUncMax->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)); } else{ Printf("This option for the calculation of the average has not been implemented yet"); return; } } } }