// Original Author: Loic Quertenmont
#include <string>
#include <vector>
#include "TROOT.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TChain.h"
#include "TObject.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TLegend.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TTree.h"
#include "TF1.h"
#include "TGraphErrors.h"
#include "TGraphAsymmErrors.h"
#include "TMultiGraph.h"
#include "TPaveText.h"
#include "TRandom3.h"
#include "tdrstyle.C"
#include "Analysis_CommonFunction.h"
#include "Analysis_Global.h"
#include "Analysis_PlotFunction.h"
#include "Analysis_PlotStructure.h"
#include "Analysis_Samples.h"
using namespace std;
/////////////////////////// CODE PARAMETERS /////////////////////////////
void Analysis_Step4(std::string InputPattern)
{
if(InputPattern=="COMPILE")return;
setTDRStyle();
gStyle->SetPadTopMargin (0.06);
gStyle->SetPadBottomMargin(0.12);
gStyle->SetPadRightMargin (0.16);
gStyle->SetPadLeftMargin (0.14);
gStyle->SetTitleSize(0.04, "XYZ");
gStyle->SetTitleXOffset(1.1);
gStyle->SetTitleYOffset(1.45);
gStyle->SetPalette(1);
gStyle->SetNdivisions(505);
unsigned int NPseudoExp = 100; //Number of PseudoExperiment to run
//7TeV DXY/DZ/ANGLE 85/326 86/327 10/251
double CosmicVetoInEfficiency7TeV = 0.26 * 0.26 * 0.04 ;
double CosmicVetoInEfficiency7TeVErr = sqrt(pow(0.03*0.26*0.04,2) + pow(0.26*0.03*0.04,2) +pow(0.26*0.04*0.01,2) + pow(0.50*CosmicVetoInEfficiency7TeV,2)); //add 50% syst uncertainty
//8TeV DXY/DZ/ANGLE 19/22 0/3 0/3
double CosmicVetoInEfficiency8TeV = 0.86 * 0.26 * 0.04 ;
double CosmicVetoInEfficiency8TeVErr = sqrt(pow(0.20*0.26*0.04,2) + pow(0.86*0.03*0.04,2) + pow(0.86*0.04*0.01,2) + pow(0.50*CosmicVetoInEfficiency8TeV,2)); //add 50% syst uncertainty
double CosmicVetoInEfficiency = CosmicVetoInEfficiency7TeV;
double CosmicVetoInEfficiencyErr = CosmicVetoInEfficiency7TeVErr;
string Input = InputPattern + "Histos.root";
TFile* InputFile = new TFile(Input.c_str(), "UPDATE");
TypeMode = TypeFromPattern(InputPattern);
//Do two loops, one for the actual background prediction and one for the
//region with TOF<1
for(unsigned int S=0; S<2; S++) {
string Suffix="";
if(S==1) Suffix="_Flip";
TH1D* HCuts_Pt = (TH1D*)GetObjectFromPath(InputFile, ("HCuts_Pt" + Suffix).c_str());
TH1D* HCuts_I = (TH1D*)GetObjectFromPath(InputFile, ("HCuts_I" + Suffix).c_str());
TH1D* HCuts_TOF = (TH1D*)GetObjectFromPath(InputFile, ("HCuts_TOF" + Suffix).c_str());
TList* list = InputFile->GetListOfKeys();
for(int d=0;d<list->GetEntries();d++){
if(!list->At(d)->IsFolder())continue;
string DirName;
DirName = DirName + list->At(d)->GetName();
if(DirName.find("Cosmic")!=string::npos) continue;
if(DirName.find("7TeV")!=string::npos){
CosmicVetoInEfficiency = CosmicVetoInEfficiency7TeV;
CosmicVetoInEfficiencyErr = CosmicVetoInEfficiency7TeVErr;
}else{
CosmicVetoInEfficiency = CosmicVetoInEfficiency8TeV;
CosmicVetoInEfficiencyErr = CosmicVetoInEfficiency8TeVErr;
}
TDirectory* directory = InputFile->GetDirectory(list->At(d)->GetName());
directory->cd();
TH1D* H_A = (TH1D*)GetObjectFromPath(directory, ("H_A" + Suffix).c_str()); if(!H_A)continue; //ABCD INFO NOT SAVED IN THIS DIRECTORY --> Skip it
TH1D* H_B = (TH1D*)GetObjectFromPath(directory, ("H_B" + Suffix).c_str());
TH1D* H_C = (TH1D*)GetObjectFromPath(directory, ("H_C" + Suffix).c_str());
TH1D* H_D = (TH1D*)GetObjectFromPath(directory, ("H_D" + Suffix).c_str());
TH1D* H_E = (TH1D*)GetObjectFromPath(directory, ("H_E" + Suffix).c_str());
TH1D* H_F = (TH1D*)GetObjectFromPath(directory, ("H_F" + Suffix).c_str());
TH1D* H_G = (TH1D*)GetObjectFromPath(directory, ("H_G" + Suffix).c_str());
TH1D* H_H = (TH1D*)GetObjectFromPath(directory, ("H_H" + Suffix).c_str());
TH1D* H_A_Cosmic = (S==0) ? (TH1D*)GetObjectFromPath(directory, ("H_A_Flip" + Suffix).c_str()) : NULL;
TH1D* H_B_Cosmic = (S==0) ? (TH1D*)GetObjectFromPath(directory, ("H_B_Flip" + Suffix).c_str()) : NULL;
TH1D* H_C_Cosmic = (S==0) ? (TH1D*)GetObjectFromPath(directory, ("H_C_Flip" + Suffix).c_str()) : NULL;
TH1D* H_D_Cosmic = (S==0) ? (TH1D*)GetObjectFromPath(directory, ("H_D_Flip" + Suffix).c_str()) : NULL;
TH1D* H_B_Binned[MaxPredBins];
TH1D* H_F_Binned[MaxPredBins];
TH1D* H_H_Binned[MaxPredBins];
TH1D* H_P_Binned[MaxPredBins];
if(TypeMode==3) PredBins=6;
for(int i=0; i<PredBins; i++) {
string Version=Suffix;
char Bin[1024];
sprintf(Bin,"_%i", i);
Version.append(Bin);
H_B_Binned[i] = (TH1D*)GetObjectFromPath(directory, ("H_B_Binned" + Version).c_str());
H_F_Binned[i] = (TH1D*)GetObjectFromPath(directory, ("H_F_Binned" + Version).c_str());
H_H_Binned[i] = (TH1D*)GetObjectFromPath(directory, ("H_H_Binned" + Version).c_str());
}
TH3D* Pred_EtaP = (TH3D*)GetObjectFromPath(directory, ("Pred_EtaP" + Suffix).c_str());
TH2D* Pred_I = (TH2D*)GetObjectFromPath(directory, ("Pred_I" + Suffix).c_str());
TH2D* Pred_TOF = (TH2D*)GetObjectFromPath(directory, ("Pred_TOF" + Suffix).c_str());
TH2D* Pred_EtaB = (TH2D*)GetObjectFromPath(directory, ("Pred_EtaB" + Suffix).c_str());
TH2D* Pred_EtaS = (TH2D*)GetObjectFromPath(directory, ("Pred_EtaS" + Suffix).c_str());
TH2D* Pred_EtaS2 = (TH2D*)GetObjectFromPath(directory, ("Pred_EtaS2" + Suffix).c_str());
TH2D* H_D_DzSidebands= (TH2D*)GetObjectFromPath(directory, "H_D_DzSidebands");
TH1D* H_D_CosmicMO=NULL;
TH2D* H_D_DzSidebands_Cosmic=NULL;
TH1D* H_B_Cosmic_Binned[MaxPredBins];
TH1D* H_F_Cosmic_Binned[MaxPredBins];
TH1D* H_H_Cosmic_Binned[MaxPredBins];
if(TypeMode==3 && DirName.find("Data")!=string::npos) {
//Only 2012 sample has pure cosmic sample, as only ratio used can use 2012 sample to make 2011 cosmic prediction
string CosmicDir = "Cosmic8TeV";
//string CosmicDir = DirName.replace(0, 4, "Cosmic");
H_D_DzSidebands_Cosmic = (TH2D*)GetObjectFromPath(InputFile, (CosmicDir + "/H_D_DzSidebands").c_str());
H_D_CosmicMO = (TH1D*)GetObjectFromPath(InputFile, (CosmicDir + "/H_D" + Suffix).c_str());
for(int i=0; i<PredBins; i++) {
string Version=Suffix;
char Bin[1024];
sprintf(Bin,"_%i", i);
Version.append(Bin);
H_B_Cosmic_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, (CosmicDir + "/H_B_Binned" + Version).c_str());
H_F_Cosmic_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, (CosmicDir + "/H_F_Binned" + Version).c_str());
H_H_Cosmic_Binned[i] = (TH1D*)GetObjectFromPath(InputFile, (CosmicDir + "/H_H_Binned" + Version).c_str());
}
}
//erase histogram created at previous iteration
//directory->Delete(("Pred_P" + Suffix + ";*").c_str());
//directory->Delete(("Pred_Mass" + Suffix + ";*").c_str());
//directory->Delete(("Pred_MassTOF" + Suffix + ";*").c_str());
//directory->Delete(("Pred_MassComb" + Suffix + ";*").c_str());
//directory->Delete(("H_P" + Suffix + ";*").c_str());
//directory->Delete(("H_P_Coll" + Suffix + ";*").c_str());
//directory->Delete(("H_P_Cosmic" + Suffix + ";*").c_str());
//take data histogram to save the resulting momentum distribution
TH1D* H_P = (TH1D*)GetObjectFromPath(directory, ("H_D" + Suffix).c_str())->Clone(("H_P" + Suffix).c_str()); H_P->Reset();
TH1D* H_P_Coll = (TH1D*)H_P->Clone(("H_P_Coll" + Suffix).c_str());
TH1D* H_P_Cosmic = (TH1D*)H_P->Clone(("H_P_Cosmic" + Suffix).c_str());
TH2D* Pred_Mass = (TH2D*)GetObjectFromPath(directory, ("Mass" + Suffix).c_str())->Clone(("Pred_Mass" + Suffix).c_str()); Pred_Mass->Reset();
TH2D* Pred_MassTOF = (TH2D*)GetObjectFromPath(directory, ("MassTOF" + Suffix).c_str())->Clone(("Pred_MassTOF" + Suffix).c_str()); Pred_MassTOF->Reset();
TH2D* Pred_MassComb = (TH2D*)GetObjectFromPath(directory, ("MassComb" + Suffix).c_str())->Clone(("Pred_MassComb" + Suffix).c_str()); Pred_MassComb->Reset();
TH2D* Pred_P = (TH2D*)GetObjectFromPath(directory, ("RegionD_P" + Suffix).c_str())->Clone(("Pred_P" + Suffix).c_str()); Pred_P->Reset();
for(int i=0; i<PredBins; i++) {
string Version=Suffix;
char Bin[1024];
sprintf(Bin,"_%i", i);
Version.append(Bin);
H_P_Binned[i] = (TH1D*)H_P->Clone(("H_P_Binned" + Version).c_str());
}
printf("Making prediction for %s\n",directory->GetName());
////////////////////////////////////////////////// MAKING THE PREDICTION
for(unsigned int CutIndex=0;CutIndex<(unsigned int)HCuts_Pt->GetXaxis()->GetNbins();CutIndex++){
//if(CutIndex<86 || CutIndex>87)continue;
double A=H_A->GetBinContent(CutIndex+1); double AErr = sqrt(A);
double B=H_B->GetBinContent(CutIndex+1); double BErr = sqrt(B);
double C=H_C->GetBinContent(CutIndex+1); double CErr = sqrt(C);
double D=H_D->GetBinContent(CutIndex+1); double DErr = sqrt(D);
double E=H_E->GetBinContent(CutIndex+1); // double EErr = sqrt(E);
double F=H_F->GetBinContent(CutIndex+1); // double FErr = sqrt(F);
double G=H_G->GetBinContent(CutIndex+1); // double GErr = sqrt(G);
double H=H_H->GetBinContent(CutIndex+1); // double HErr = sqrt(H);
double A_Cosmic=1, B_Cosmic=0, C_Cosmic=0, D_Cosmic=0;
double AErr_Cosmic=0, BErr_Cosmic=0, CErr_Cosmic=0, DErr_Cosmic=0;
if(S==0 && TypeMode==5){
A_Cosmic=H_A_Cosmic->GetBinContent(CutIndex+1) * CosmicVetoInEfficiency;
B_Cosmic=H_B_Cosmic->GetBinContent(CutIndex+1) * CosmicVetoInEfficiency;
C_Cosmic=H_C_Cosmic->GetBinContent(CutIndex+1) * CosmicVetoInEfficiency;
D_Cosmic=H_D_Cosmic->GetBinContent(CutIndex+1) * CosmicVetoInEfficiency;
AErr_Cosmic=sqrt( pow(sqrt(H_A_Cosmic->GetBinContent(CutIndex+1)) * CosmicVetoInEfficiency,2) + pow(H_A_Cosmic->GetBinContent(CutIndex+1) * CosmicVetoInEfficiencyErr,2) );
BErr_Cosmic=sqrt( pow(sqrt(H_B_Cosmic->GetBinContent(CutIndex+1)) * CosmicVetoInEfficiency,2) + pow(H_B_Cosmic->GetBinContent(CutIndex+1) * CosmicVetoInEfficiencyErr,2) );
CErr_Cosmic=sqrt( pow(sqrt(H_C_Cosmic->GetBinContent(CutIndex+1)) * CosmicVetoInEfficiency,2) + pow(H_C_Cosmic->GetBinContent(CutIndex+1) * CosmicVetoInEfficiencyErr,2) );
DErr_Cosmic=sqrt( pow(sqrt(H_D_Cosmic->GetBinContent(CutIndex+1)) * CosmicVetoInEfficiency,2) + pow(H_D_Cosmic->GetBinContent(CutIndex+1) * CosmicVetoInEfficiencyErr,2) );
if(CutIndex==44){
printf("scale factor = %f+-%f\n", CosmicVetoInEfficiency, CosmicVetoInEfficiencyErr);
printf("%E+-%E %E+-%E %E+-%E\n", A_Cosmic, AErr_Cosmic, B_Cosmic, BErr_Cosmic, C_Cosmic, CErr_Cosmic);
}
A = A - A_Cosmic; AErr = sqrt(AErr*AErr + AErr_Cosmic*AErr_Cosmic);
B = B - B_Cosmic; BErr = sqrt(BErr*BErr + BErr_Cosmic*BErr_Cosmic);
C = C - C_Cosmic; CErr = sqrt(CErr*CErr + CErr_Cosmic*CErr_Cosmic);
//D = D - D_Cosmic; DErr = sqrt(DErr*DErr + DErr_Cosmic*DErr_Cosmic);
}
double B_Binned[MaxPredBins];
double F_Binned[MaxPredBins];
double H_Binned[MaxPredBins];
//double P_Binned[MaxPredBins];
//double Perr_Binned[MaxPredBins];
double B_Cosmic_Binned[MaxPredBins];
double F_Cosmic_Binned[MaxPredBins];
double H_Cosmic_Binned[MaxPredBins];
for(int i=0; i<PredBins; i++) {
B_Binned[i]=H_B_Binned[i]->GetBinContent(CutIndex+1);
F_Binned[i]=H_F_Binned[i]->GetBinContent(CutIndex+1);
H_Binned[i]=H_H_Binned[i]->GetBinContent(CutIndex+1);
//P_Binned[i]=0;
//Perr_Binned[i]=0;
B_Cosmic_Binned[i]=H_B_Cosmic_Binned[i]->GetBinContent(CutIndex+1);
F_Cosmic_Binned[i]=H_F_Cosmic_Binned[i]->GetBinContent(CutIndex+1);
H_Cosmic_Binned[i]=H_H_Cosmic_Binned[i]->GetBinContent(CutIndex+1);
}
double P=0;
double Perr=0;
double P_Coll=0;
double Perr_Coll=0;
double P_Cosmic=0;
double Perr_Cosmic=0;
printf("%4i --> Pt>%7.2f I>%6.2f TOF>%+5.2f --> A=%6.2E B=%6.2E C=%6.2E D=%6.2E E=%6.2E F=%6.2E G=%6.2E H=%6.2E",CutIndex,HCuts_Pt->GetBinContent(CutIndex+1), HCuts_I->GetBinContent(CutIndex+1), HCuts_TOF->GetBinContent(CutIndex+1),A, B, C, D, E, F, G, H );
if(E>0){
//Prediction in Pt-Is-TOF plane
P = (A*F*G)/(E*E);
Perr = sqrt( ((pow(F*G,2)* A + pow(A*G,2)*F + pow(A*F,2)*G)/pow(E,4)) + (pow((2*A*F*G)/pow(E,3),2)*E));
}else if(A>0){
//Prediction in Pt-Is plane
P = ((C*B)/A);
Perr = sqrt( pow(CErr*B/A,2) + pow(BErr*C/A,2) + pow((AErr*B*C/(A*A)),2) );
if(S==0 && TypeMode==5){
P_Coll = P;
Perr_Coll = Perr;
P_Cosmic = -1;
if(A_Cosmic>0 && B_Cosmic>0 && C_Cosmic>0){
P_Cosmic = ((C_Cosmic*B_Cosmic)/A_Cosmic);
Perr_Cosmic = sqrt( pow(CErr_Cosmic*B_Cosmic/A_Cosmic,2) + pow(BErr_Cosmic*C_Cosmic/A_Cosmic,2) + pow(AErr_Cosmic*B_Cosmic*C_Cosmic/(A_Cosmic*A_Cosmic),2) );
}else if(D_Cosmic>0){
P_Cosmic = D_Cosmic * CosmicVetoInEfficiency;
Perr_Cosmic = sqrt( pow(DErr_Cosmic * CosmicVetoInEfficiency,2) + pow(D_Cosmic*CosmicVetoInEfficiencyErr,2));
}else{
P_Cosmic = 3/2.0 * CosmicVetoInEfficiency;
Perr_Cosmic = P_Cosmic; //100% uncertainty
}
P = P_Coll + P_Cosmic;
Perr = sqrt( Perr_Coll*Perr_Coll + Perr_Cosmic*Perr_Cosmic);
}
}else if(F>0){
//Predict the number of cosmics passing all cuts as number passing in dz sideband times the ratio of tracks in the sideband
//vs number in central region as determined by pure cosmic sample
//Multile sidebands are made to check for background consistency, the fifth one is used for the actual prediction
double D_Sideband = 0;
double D_Sideband_Cosmic = 0;
if(DirName.find("Data")!=string::npos) {
D_Sideband = H_D_DzSidebands->GetBinContent(CutIndex+1, 5);
double D_Cosmic = H_D_CosmicMO->GetBinContent(CutIndex+1);
D_Sideband_Cosmic = H_D_DzSidebands_Cosmic->GetBinContent(CutIndex+1, 5);
if(D_Sideband_Cosmic>0) {
P_Cosmic = D_Sideband * D_Cosmic / D_Sideband_Cosmic;
Perr_Cosmic = sqrt( (pow(D_Cosmic/D_Sideband_Cosmic,2)*D_Sideband) + (pow(D_Sideband/D_Sideband_Cosmic,2)*D_Cosmic) + (pow((D_Cosmic*(D_Sideband)/(D_Sideband_Cosmic*D_Sideband_Cosmic)),2)*D_Sideband_Cosmic) );
}
}
//Prediction in Pt-TOF plane
for(int i=0; i<PredBins; i++) {
//Subtract the expected cosmic tracks from each region
double B_Bin = B_Binned[i] - B_Cosmic_Binned[i]*D_Sideband/D_Sideband_Cosmic;
double F_Bin = F_Binned[i] - F_Cosmic_Binned[i]*D_Sideband/D_Sideband_Cosmic;
double H_Bin = H_Binned[i] - H_Cosmic_Binned[i]*D_Sideband/D_Sideband_Cosmic;
double Berr = sqrt(B_Binned[i] + (pow(B_Cosmic_Binned[i]/D_Sideband_Cosmic,2)*D_Sideband) + (pow(D_Sideband/D_Sideband_Cosmic,2)*B_Cosmic_Binned[i]) + (pow((B_Cosmic_Binned[i]*(D_Sideband)/(D_Sideband_Cosmic*D_Sideband_Cosmic)),2)*D_Sideband_Cosmic) );
double Ferr = sqrt(F_Binned[i] + (pow(F_Cosmic_Binned[i]/D_Sideband_Cosmic,2)*D_Sideband) + (pow(D_Sideband/D_Sideband_Cosmic,2)*F_Cosmic_Binned[i]) + (pow((F_Cosmic_Binned[i]*(D_Sideband)/(D_Sideband_Cosmic*D_Sideband_Cosmic)),2)*D_Sideband_Cosmic) );
double Herr = sqrt(H_Binned[i] + (pow(H_Cosmic_Binned[i]/D_Sideband_Cosmic,2)*D_Sideband) + (pow(D_Sideband/D_Sideband_Cosmic,2)*H_Cosmic_Binned[i]) + (pow((H_Cosmic_Binned[i]*(D_Sideband)/(D_Sideband_Cosmic*D_Sideband_Cosmic)),2)*D_Sideband_Cosmic) );
double P_Binned = ((H_Bin*B_Bin)/F_Bin);
double Perr_Binned = (pow(Berr/Ferr,2)*Herr) + (pow(Herr/Ferr,2)*Berr) + (pow((Berr*(Herr)/(Ferr*Ferr)),2)*Ferr);
H_P_Binned[i]->SetBinContent(CutIndex+1, P_Binned);
H_P_Binned[i]->SetBinError(CutIndex+1, sqrt(Perr_Binned));
P_Coll += P_Binned;
Perr_Coll += Perr_Binned;
}
Perr_Coll = sqrt(Perr_Coll);
P = P_Coll + P_Cosmic;
Perr = sqrt(Perr_Coll*Perr_Coll + Perr_Cosmic*Perr_Cosmic);
//Add in systematic contribution, now done in step 6
//Perr = sqrt(Perr*Perr + P_Coll*P_Coll*0.2*0.2 + P_Cosmic*P_Cosmic*0.8*0.8);
}else if(G>0){
//Prediction in Ias-TOF plane
P = ((C*H)/G);
Perr = sqrt( (pow(H/G,2)*C) + (pow(C/G,2)*H) + (pow((H*(C)/(G*G)),2)*G) );
}
H_P_Coll->SetBinContent(CutIndex+1, P_Coll);
H_P_Coll->SetBinError (CutIndex+1, Perr_Coll);
H_P_Cosmic->SetBinContent(CutIndex+1,P_Cosmic);
H_P_Cosmic->SetBinError (CutIndex+1,Perr_Cosmic);
H_P->SetBinContent(CutIndex+1,P);
H_P->SetBinError (CutIndex+1,Perr);
if(P==0 || isnan((float)P)) {printf("\n"); continue;} //Skip this CutIndex --> No Prediction possible
printf(" --> D=%6.2E vs Pred = %6.2E +- %6.2E (%6.2E%%)\n", D, P, Perr, 100.0*Perr/P );
if(TypeMode>2)continue; //Need to compute mass predicted distribution ONLY for TkOnly and TkTOF
TH1D* Pred_EtaB_Proj = Pred_EtaB ->ProjectionY("ProjEtaB" ,CutIndex+1,CutIndex+1);
TH1D* Pred_EtaS_Proj = Pred_EtaS ->ProjectionY("ProjEtaS" ,CutIndex+1,CutIndex+1);
TH1D* Pred_EtaS2_Proj = Pred_EtaS2->ProjectionY("ProjEtaS2",CutIndex+1,CutIndex+1);
TH1D* Pred_EtaB_Proj_PE = (TH1D*)Pred_EtaB_Proj ->Clone("Pred_EtaB_Proj_PE"); Pred_EtaB_Proj_PE ->Reset();
TH1D* Pred_EtaS_Proj_PE = (TH1D*)Pred_EtaS_Proj ->Clone("Pred_EtaS_Proj_PE"); Pred_EtaS_Proj_PE ->Reset();
TH1D* Pred_EtaS2_Proj_PE = (TH1D*)Pred_EtaS2_Proj->Clone("Pred_EtaS2_Proj_PE"); Pred_EtaS2_Proj_PE->Reset();
Pred_EtaP->GetXaxis()->SetRange(CutIndex+1,CutIndex+1);
TH2D* Pred_EtaPWeighted = (TH2D*)Pred_EtaP->Project3D("zy");
TH2D* Pred_EtaPWeighted_PE = (TH2D*)Pred_EtaPWeighted->Clone("Pred_EtaPWeightedPE"); Pred_EtaPWeighted_PE->Reset();
TH1D* Pred_I_Proj = Pred_I->ProjectionY("ProjI",CutIndex+1,CutIndex+1);
TH1D* Pred_T_Proj = Pred_TOF->ProjectionY("ProjT",CutIndex+1,CutIndex+1);
TH1D* Pred_I_ProjPE = (TH1D*) Pred_I_Proj->Clone("Pred_I_ProjPE"); Pred_I_ProjPE->Reset();
TH1D* Pred_T_ProjPE = (TH1D*) Pred_T_Proj->Clone("Pred_T_ProjPE"); Pred_T_ProjPE->Reset();
TH2D* Pred_Prof_Mass = new TH2D("Pred_Prof_Mass" ,"Pred_Prof_Mass" ,MassNBins,0,MassHistoUpperBound, NPseudoExp, 0, NPseudoExp);
TH2D* Pred_Prof_MassTOF = new TH2D("Pred_Prof_MassTOF" ,"Pred_Prof_MassTOF" ,MassNBins,0,MassHistoUpperBound, NPseudoExp, 0, NPseudoExp);
TH2D* Pred_Prof_MassComb = new TH2D("Pred_Prof_MassComb","Pred_Prof_MassComb",MassNBins,0,MassHistoUpperBound, NPseudoExp, 0, NPseudoExp);
for(int x=0;x<Pred_Mass->GetNbinsY()+1;x++){
for(unsigned int pe=0;pe<NPseudoExp;pe++){
Pred_Prof_Mass ->SetBinContent(x, pe, 0);
Pred_Prof_MassTOF ->SetBinContent(x, pe, 0);
Pred_Prof_MassComb->SetBinContent(x, pe, 0);
}
}
TRandom3* RNG = new TRandom3();
printf("Predicting (%4i / %4i) :",CutIndex+1,HCuts_Pt->GetXaxis()->GetNbins());
int TreeStep = NPseudoExp/50;if(TreeStep==0)TreeStep=1;
for(unsigned int pe=0;pe<NPseudoExp;pe++){
if(pe%TreeStep==0){printf(".");fflush(stdout);}
TH1D* tmpH_Mass = new TH1D("tmpH_Mass" ,"tmpH_Mass" ,MassNBins,0,MassHistoUpperBound);
TH1D* tmpH_MassTOF = new TH1D("tmpH_MassTOF" ,"tmpH_MassTOF" ,MassNBins,0,MassHistoUpperBound);
TH1D* tmpH_MassComb = new TH1D("tmpH_MassComb","tmpH_MassComb",MassNBins,0,MassHistoUpperBound);
double PE_A=RNG->Poisson(A);
double PE_B=RNG->Poisson(B);
double PE_C=RNG->Poisson(C);
double PE_E=RNG->Poisson(E);
double PE_F=RNG->Poisson(F);
double PE_G=RNG->Poisson(G);
double PE_P = 0;
if(E>0){ PE_P = (PE_E>0 ? (PE_A*PE_F*PE_G)/(PE_E*PE_E) : 0);
}else if(A>0){ PE_P = (PE_A>0 ? ((PE_C*PE_B)/PE_A) : 0);
}
for(int i=0;i<Pred_EtaB_Proj_PE->GetNbinsX()+1;i++){Pred_EtaB_Proj_PE->SetBinContent(i,RNG->Poisson(Pred_EtaB_Proj->GetBinContent(i)) );} Pred_EtaB_Proj_PE->Scale(1.0/Pred_EtaB_Proj_PE->Integral());
for(int i=0;i<Pred_EtaS_Proj_PE->GetNbinsX()+1;i++){Pred_EtaS_Proj_PE->SetBinContent(i,RNG->Poisson(Pred_EtaS_Proj->GetBinContent(i)) );} Pred_EtaS_Proj_PE->Scale(1.0/Pred_EtaS_Proj_PE->Integral());
for(int i=0;i<Pred_EtaS2_Proj_PE->GetNbinsX()+1;i++){Pred_EtaS2_Proj_PE->SetBinContent(i,RNG->Poisson(Pred_EtaS2_Proj->GetBinContent(i)) );} Pred_EtaS2_Proj_PE->Scale(1.0/Pred_EtaS2_Proj_PE->Integral());
for(int i=0;i<Pred_EtaPWeighted_PE->GetNbinsX()+1;i++){
for(int j=0;j<Pred_EtaPWeighted_PE->GetNbinsY()+1;j++){
Pred_EtaPWeighted_PE->SetBinContent(i,j,RNG->Poisson(Pred_EtaPWeighted->GetBinContent(i,j)));
}}
double WeightP = 0.0;
for(int x=0;x<=Pred_EtaPWeighted_PE->GetXaxis()->GetNbins();x++){
WeightP = 0.0;
if(Pred_EtaB_Proj_PE->GetBinContent(x)>0){
WeightP = Pred_EtaS_Proj_PE ->GetBinContent(x)/Pred_EtaB_Proj_PE->GetBinContent(x);
if(TypeMode==2)WeightP*= Pred_EtaS2_Proj_PE->GetBinContent(x)/Pred_EtaB_Proj_PE->GetBinContent(x);
}
for(int y=0;y<=Pred_EtaPWeighted_PE->GetYaxis()->GetNbins();y++){
Pred_EtaPWeighted_PE->SetBinContent(x,y,Pred_EtaPWeighted_PE->GetBinContent(x,y)*WeightP);
}
}
TH1D* Pred_P_ProjPE = Pred_EtaPWeighted_PE->ProjectionY("Pred_P_ProjPE"); Pred_P_ProjPE->Scale(1.0/Pred_P_ProjPE->Integral());
for(int i=0;i<Pred_I_ProjPE->GetNbinsX()+1;i++){Pred_I_ProjPE->SetBinContent(i,RNG->Poisson(Pred_I_Proj->GetBinContent(i)) );} Pred_I_ProjPE->Scale(1.0/Pred_I_ProjPE->Integral());
for(int i=0;i<Pred_T_ProjPE->GetNbinsX()+1;i++){Pred_T_ProjPE->SetBinContent(i,RNG->Poisson(Pred_T_Proj->GetBinContent(i)) );} Pred_T_ProjPE->Scale(1.0/Pred_T_ProjPE->Integral());
//save the predP distribution
for(int x=0;x<Pred_P_ProjPE->GetNbinsX()+1;x++){Pred_P->SetBinContent(CutIndex+1, x, Pred_P->GetBinContent(CutIndex+1, x) + Pred_P_ProjPE->GetBinContent(x) * PE_P);};
double Proba, MI, MComb;//, MT=0, ProbaT=0;
for(int x=0;x<Pred_P_ProjPE->GetNbinsX()+1;x++){ if(Pred_P_ProjPE->GetBinContent(x)<=0.0){continue;} const double& p = Pred_P_ProjPE->GetBinCenter(x);
for(int y=0;y<Pred_I_ProjPE->GetNbinsX()+1;y++){ if(Pred_I_ProjPE->GetBinContent(y)<=0.0){continue;} const double& i = Pred_I_ProjPE->GetBinCenter(y);
Proba = Pred_P_ProjPE->GetBinContent(x) * Pred_I_ProjPE->GetBinContent(y); if(Proba<=0 || isnan((float)Proba))continue;
MI = GetMass(p,i, false);
MComb = MI;
tmpH_Mass->Fill(MI,Proba);
//if(MI>500 && MI<800 && pe==0 && (CutIndex==3 || CutIndex==4)){printf("Index=%i p=%7.2f i=%7.2f M=%7.2f P=%7.5E( = %7.5E x %7.5E) --> %7.5E | %7.5E\n",CutIndex,p,i,MI,Proba,Pred_P_ProjPE->GetBinContent(x),Pred_I_ProjPE->GetBinContent(y),Pred_P_ProjPE->GetBinContent(x)*1,Pred_I_ProjPE->GetBinContent(y)*1);}
//commented part there is related to the prediction of the mass reconstructed from TOF.
//if(TypeMode==2){
//for(int z=0;z<Pred_T_ProjPE->GetNbinsX()+1;z++){ if(Pred_T_ProjPE->GetBinContent(z)<=0.0){continue;} const double& t = Pred_T_ProjPE->GetBinCenter(z);
// ProbaT = Proba * Pred_T_ProjPE->GetBinContent(z); if(ProbaT<=0 || isnan(ProbaT))continue;
// MT = GetTOFMass(p,t);
// tmpH_MassTOF->Fill(MT,ProbaT);
// MComb = GetMassFromBeta(p, (GetIBeta(i, false) + (1/t))*0.5 );
// tmpH_MassComb->Fill(MComb,ProbaT);
//}}else{
tmpH_MassComb->Fill(MComb,Proba);
//}
}}
for(int x=0;x<tmpH_Mass->GetNbinsX()+1;x++){
Pred_Prof_Mass ->SetBinContent(x, pe, tmpH_Mass ->GetBinContent(x) * PE_P);
Pred_Prof_MassTOF ->SetBinContent(x, pe, tmpH_MassTOF ->GetBinContent(x) * PE_P);
Pred_Prof_MassComb->SetBinContent(x, pe, tmpH_MassComb->GetBinContent(x) * PE_P);
if(isnan((float)(tmpH_Mass ->GetBinContent(x) * PE_P))){printf("%f x %f\n",tmpH_Mass ->GetBinContent(x),PE_P); fflush(stdout);exit(0);}
}
delete Pred_P_ProjPE;
delete tmpH_Mass;
delete tmpH_MassTOF;
delete tmpH_MassComb;
}printf("\n");
for(int x=0;x<Pred_Mass->GetNbinsY()+1;x++){
double Mean=0, MeanTOF=0, MeanComb=0;
for(unsigned int pe=0;pe<NPseudoExp;pe++){
Mean += Pred_Prof_Mass ->GetBinContent(x, pe);
MeanTOF += Pred_Prof_MassTOF ->GetBinContent(x, pe);
MeanComb += Pred_Prof_MassComb->GetBinContent(x, pe);
}Mean/=NPseudoExp; MeanTOF/=NPseudoExp; MeanComb/=NPseudoExp;
double Err=0, ErrTOF=0, ErrComb=0;
for(unsigned int pe=0;pe<NPseudoExp;pe++){
Err += pow(Mean - Pred_Prof_Mass ->GetBinContent(x, pe),2);
ErrTOF += pow(MeanTOF - Pred_Prof_MassTOF ->GetBinContent(x, pe),2);
ErrComb += pow(MeanComb - Pred_Prof_MassComb->GetBinContent(x, pe),2);
}Err=sqrt(Err/(NPseudoExp-1)); ErrTOF=sqrt(ErrTOF/(NPseudoExp-1)); ErrComb=sqrt(ErrComb/(NPseudoExp-1));
Pred_Mass ->SetBinContent(CutIndex+1,x,Mean ); Pred_Mass ->SetBinError(CutIndex+1,x,Err );
Pred_MassTOF ->SetBinContent(CutIndex+1,x,MeanTOF ); Pred_MassTOF ->SetBinError(CutIndex+1,x,ErrTOF );
Pred_MassComb->SetBinContent(CutIndex+1,x,MeanComb); Pred_MassComb ->SetBinError(CutIndex+1,x,ErrComb);
}
delete Pred_EtaB_Proj_PE;
delete Pred_EtaS_Proj_PE;
delete Pred_EtaS2_Proj_PE;
delete Pred_Prof_Mass;
delete Pred_Prof_MassTOF;
delete Pred_Prof_MassComb;
delete Pred_EtaPWeighted_PE;
delete Pred_I_ProjPE;
delete Pred_T_ProjPE;
delete Pred_I_Proj;
delete Pred_T_Proj;
delete Pred_EtaB_Proj;
delete Pred_EtaS_Proj;
delete Pred_EtaS2_Proj;
delete Pred_EtaPWeighted;
}
//scale it down by the number of PseudoExperiment to get right normalization
Pred_P->Scale(1.0/NPseudoExp);
//save histogram to file
Pred_P ->Write();
H_P ->Write();
Pred_Mass ->Write();
Pred_MassTOF ->Write();
Pred_MassComb->Write();
H_P_Coll->Write();
H_P_Cosmic->Write();
if(TypeMode==3) {
for(int i=0; i<PredBins; i++) {
H_P_Binned[i]->Write();
}
}
//directory->Delete("H_P;1");
////////////////////////////////////////////////// DUMP USEFUL INFORMATION
FILE* pFile = fopen((InputPattern+"/Info_"+directory->GetName()+Suffix+".txt").c_str(),"w");
for(unsigned int CutIndex=0;CutIndex<(unsigned int)HCuts_Pt->GetXaxis()->GetNbins();CutIndex++){
const double& A=H_A->GetBinContent(CutIndex+1);
const double& B=H_B->GetBinContent(CutIndex+1);
const double& C=H_C->GetBinContent(CutIndex+1);
const double& D=H_D->GetBinContent(CutIndex+1);
const double& E=H_E->GetBinContent(CutIndex+1);
const double& F=H_F->GetBinContent(CutIndex+1);
const double& G=H_G->GetBinContent(CutIndex+1);
const double& H=H_H->GetBinContent(CutIndex+1);
fprintf(pFile ,"CutIndex=%4i --> (Pt>%6.2f I>%6.3f TOF>%6.3f) Ndata=%+6.2E NPred=%6.3E+-%6.3E (=%6.3E+-%6.3E + %6.3E+-%6.3E) <--> A=%6.2E B=%6.2E C=%6.2E D=%6.2E E=%6.2E F=%6.2E G=%6.2E H=%6.2E\n",CutIndex,HCuts_Pt ->GetBinContent(CutIndex+1), HCuts_I ->GetBinContent(CutIndex+1), HCuts_TOF->GetBinContent(CutIndex+1), D,H_P->GetBinContent(CutIndex+1),H_P->GetBinError(CutIndex+1), H_P_Coll->GetBinContent(CutIndex+1),H_P_Coll->GetBinError(CutIndex+1), H_P_Cosmic->GetBinContent(CutIndex+1),H_P_Cosmic->GetBinError(CutIndex+1), A, B, C, D, E, F, G, H);
}
fprintf(pFile,"--------------------\n");
fclose(pFile);
}//end loop on sub directory
}//End of loop on two predictions
}