diff --git a/createDelphesScripts.sh b/createDelphesScripts.sh new file mode 100755 index 00000000..3a34b4d9 --- /dev/null +++ b/createDelphesScripts.sh @@ -0,0 +1,60 @@ +#!/bin/bash + +# Create fast simulation resolutions plotting scripts and submit to ifarm +METHOD="Ele" # Switch this to select reconstruction method from {"Ele","JB","DA"} +script="$PWD/macro/analysis_resolution.C" +postScript="$PWD/macro/postprocess_resolution.C" +submitScript="$PWD/submit.sh" +jobScript="$PWD/job.sh" +out="$PWD/macro/" +cd datarec +for file in *-xm25.config +do + energies=`echo $file | sed "s;.config;;g"` + config="${METHOD}_${energies}" + mkdir -p $out/$config + cp $script $out/$config/ + newscript=${out}${config}/*.C + eleIn=`echo $file | grep -Eo "[0-9][0-9]*x[0-9]*" | sed "s;x.*;;g"` + beamIn=`echo $file | grep -Eo "[0-9][0-9]*x[0-9]*" | sed "s;.*x;;g"` + xAng=`echo $file | grep -Eo "[0-9][0-9]*x[0-9]*-x[0-9]*" | sed "s;.*-x;;g"` + xAngM=`echo $file | grep -Eo "[0-9][0-9]*x[0-9]*-xm[0-9]*" | sed "s;.*-xm;;g"` + echo "file=${file}" + echo "energies=${energies}" + echo "config=${config}" + echo "newscript=${newscript}" + echo "eleIn=${eleIn}" + echo "beamIn=${beamIn}" + echo "xAng=${xAng}" + echo "xAngM=${xAngM}" + + # Modify analysis script + sed -i "s;datarec/dis-5x41;datarec/${energies};g" $newscript + sed -i "s;Ele_dis-5x41;${config};g" $newscript + sed -i "s;eleBeamEn=5;eleBeamEn=${eleIn};g" $newscript + sed -i "s;ionBeamEn=41;ionBeamEn=${beamIn};g" $newscript + if [ $xAng ]; then + sed -i "s;crossingAngle=25;crossingAngle=-${xAng};g" $newscript + fi + if [ $xAngM ]; then + sed -i "s;crossingAngle=25;crossingAngle=${xAngM};g" $newscript + fi + sed -i "s;Ele;${METHOD};g" $newscript + + # Postprocessor + cp $postScript $out/$config + sed -i "s;dis-5x41;${config};g" $out/$config/postprocess*.C + sed -i "s;testheader;${eleIn}x${beamIn}GeV;g" $out/$config/postprocess*.C + + # And job scripts + cp $submitScript $out/$config + cp $jobScript $out/$config + sed -i "s;dis-5x41;${config};g" $out/$config/*.sh + + # And submit + sbatch $out/$config/submit.sh + echo -------------------- + +done +cd .. +echo DONE diff --git a/createFullScripts.sh b/createFullScripts.sh new file mode 100755 index 00000000..7a6ca8ce --- /dev/null +++ b/createFullScripts.sh @@ -0,0 +1,48 @@ +#!/bin/bash +# Create full simulation resolution plotting scripts and submit to ifarm +script="$PWD/macro/analysis_resolution_SD_Full.C" +postScript="$PWD/macro/postprocess_resolution_SD.C" +submitScript="$PWD/submit.sh" +jobScript="$PWD/job.sh" +out="$PWD/macro/canyonlands/" +cd datarec/canyonlands/ +for file in */*.config +do + config=`echo $file | sed "s;.config;;g" | grep -Eo 'canyonlands/.*' | sed 's;canyonlands/;;g' | sed 's;/files;;g'` + mkdir -p $out/$config + cp $script $out/$config/ + newscript=${out}${config}/*.C + eleIn=`echo $file | grep -Eo "[0-9][0-9]*x[0-9]*" | sed "s;x.*;;g"` + beamIn=`echo $file | grep -Eo "[0-9][0-9]*x[0-9]*" | sed "s;.*x;;g"` + xAng=`echo $file | grep -Eo "[0-9][0-9]*x[0-9]*-x[0-9]*" | sed "s;.*-x;;g"` + xAngM=`echo $file | grep -Eo "[0-9][0-9]*x[0-9]*-xm[0-9]*" | sed "s;.*-xm;;g"` + echo "name=${name}" + echo "file=${file}" + echo "config=${config}" + echo "newscript=${newscript}" + echo "eleIn=${eleIn}" + echo "beamIn=${beamIn}" + echo "xAng=${xAng}" + echo "xAngM=${xAngM}" + + sed -i "s;datarec/dis-5x41.config;${file};g" $newscript + sed -i "s;full_dis-5x41;full_${config};g" $newscript + sed -i "s;eleBeamEn=5;eleBeamEn=${eleIn};g" $newscript + sed -i "s;ionBeamEn=41;ionBeamEn=${beamIn};g" $newscript + if [ $xAng ]; then + sed -i "s;crossingAngle=25;crossingAngle=-${xAng};g" $newscript + fi + if [ $xAngM ]; then + sed -i "s;icrossingAngle=25;crossingAngle=${xAngM};g" $newscript + fi + cp $postScript $out/$config + sed -i "s;dis-5x41;${config};g" $out/$config/*.C + sed -i "s;testheader;${eleIn}x${beamIn}GeV;g" $out/$config/*.C + cp $submitScript $out/$config + cp $jobScript $out/$config + sed -i "s;dis-5x41;canyonlands/${config};g" $out/$config/*.sh + sbatch $out/$config/submit.sh + echo -------------------- +done +cd .. +echo DONE diff --git a/datarec/dis-5x41-xm25.config b/datarec/dis-5x41-xm25.config new file mode 100644 index 00000000..55df4409 --- /dev/null +++ b/datarec/dis-5x41-xm25.config @@ -0,0 +1,3 @@ +1 4.002e-4 /lustre19/expphy/volatile/halla/solid/dcbyer/eic-data/new/dis-5x41-xm25.root -1 +100 2.869e-7 /lustre19/expphy/volatile/halla/solid/dcbyer/eic-data/new/dis-5x41-q2min100-xm25.root -1 + diff --git a/job.sh b/job.sh new file mode 100755 index 00000000..3d0fcd0f --- /dev/null +++ b/job.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +cd /work/clas12/users/$USER/eic3/largex-eic +echo "root -q -b macro/dis-5x41/analysis_resolution.C" | ./container/shell.sh +echo "root -q -b macro/dis-5x41/postprocess_resolution.C" | ./container/shell.sh +echo DONE diff --git a/macro/analysis_barak.C b/macro/analysis_barak.C new file mode 100644 index 00000000..7ff16c32 --- /dev/null +++ b/macro/analysis_barak.C @@ -0,0 +1,47 @@ +R__LOAD_LIBRARY(Largex) + +// make resolution plots +void analysis_barak( + TString infiles="datarec/dis-18x275-xm25.config", /* list of input files */ + Double_t eleBeamEn=18, /* electron beam energy [GeV] */ + Double_t ionBeamEn=275, /* ion beam energy [GeV] */ + Double_t crossingAngle=25, /* crossing angle [mrad] */ + TString outfilePrefix="barak_dis-18x275-xm25" /* output filename prefix*/ +) { + + // setup analysis ======================================== + AnalysisDelphes *A = new AnalysisDelphes( + infiles, + eleBeamEn, + ionBeamEn, + crossingAngle, + outfilePrefix + ); + A->NBINS = 1; // use this to set the number of bins along each axis, e.g., z binning (except resolution axes) for each overall bin in e.g. x and Q2 + A->NBINSRES = 100; // use this to set the number of bins along the resolution axes for each overall bin in e.g. x and Q2 + // A->RESHIGH = 1; + // A->RESLOW = -1; + //A->maxEvents = 30000; // use this to limit the number of events + A->SetReconMethod("Ele"); // set reconstruction method + // A->SetReconMethod("JB"); // set reconstruction method + //A->SetReconMethod("DA"); // set reconstruction method + A->AddFinalState("pipTrack"); // pi+ final state + A->AddFinalState("pimTrack"); // pi- final state + A->AddFinalState("KpTrack"); // K+ final state + A->AddFinalState("KmTrack"); // K- final state + A->AddFinalState("pTrack"); // proton final state + + // define cuts ==================================== + A->AddBinScheme("w"); A->BinScheme("w")->BuildBin("Min",3.0); // W > 3 GeV + A->AddBinScheme("y"); A->BinScheme("y")->BuildBin("Range",0.01,0.95); // 0.01 < y < 0.95 + A->AddBinScheme("z"); A->BinScheme("z")->BuildBin("Range",0.2,0.9); // 0.2 < z < 0.9 + A->AddBinScheme("xF"); A->BinScheme("xF")->BuildBin("Min",0.0); // xF > 0 + A->AddBinScheme("ptLab"); A->BinScheme("ptLab")->BuildBin("Min",0.1); // pT_lab > 0.1 GeV (tracking limit) + + // set binning scheme ==================================== + A->AddBinScheme("q2"); A->BinScheme("q2")->BuildBins( 25, 0.1, 10000, true ); + A->AddBinScheme("x"); A->BinScheme("x")->BuildBins( 25, 0.00001, 1, true ); + + // perform the analysis ================================== + A->Execute(); +}; diff --git a/macro/analysis_resolution.C b/macro/analysis_resolution.C new file mode 100644 index 00000000..aa408eb1 --- /dev/null +++ b/macro/analysis_resolution.C @@ -0,0 +1,45 @@ +R__LOAD_LIBRARY(Largex) + +// make resolution plots +void analysis_resolution( + TString infiles="datarec/dis-5x41.config", /* list of input files */ + Double_t eleBeamEn=5, /* electron beam energy [GeV] */ + Double_t ionBeamEn=41, /* ion beam energy [GeV] */ + Double_t crossingAngle=25, /* crossing angle [mrad] */ + TString outfilePrefix="Ele_dis-5x41" /* output filename prefix*/ +) { + + // setup analysis ======================================== + AnalysisDelphes *A = new AnalysisDelphes( + infiles, + eleBeamEn, + ionBeamEn, + crossingAngle, + outfilePrefix + ); + A->NBINS = 1; // use this to set the number of bins along each axis, e.g., z binning (except resolution axes) for each overall bin in e.g. x and Q2 + A->NBINSRES = 100; // use this to set the number of bins along the resolution axes for each overall bin in e.g. x and Q2 + // A->RESHIGH = 1; + // A->RESLOW = -1; + //A->maxEvents = 30000; // use this to limit the number of events + A->SetReconMethod("Ele"); // set reconstruction method + A->AddFinalState("pipTrack"); // pion final state + // A->AddFinalState("pimTrack"); // pion final state + // A->AddFinalState("KpTrack"); // pion final state + // A->AddFinalState("KmTrack"); // pion final state + // A->AddFinalState("pTrack"); // pion final state + + // define cuts ==================================== + A->AddBinScheme("w"); A->BinScheme("w")->BuildBin("Min",3.0); // W > 3 GeV + A->AddBinScheme("y"); A->BinScheme("y")->BuildBin("Range",0.01,0.95); // 0.01 < y < 0.95 + A->AddBinScheme("z"); A->BinScheme("z")->BuildBin("Range",0.2,0.9); // 0.2 < z < 0.9 + A->AddBinScheme("xF"); A->BinScheme("xF")->BuildBin("Min",0.0); // xF > 0 + A->AddBinScheme("ptLab"); A->BinScheme("ptLab")->BuildBin("Min",0.1); // pT_lab > 0.1 GeV (tracking limit) + + // set binning scheme ==================================== + A->AddBinScheme("q2"); A->BinScheme("q2")->BuildBins( 4, 1, 1000, true ); + A->AddBinScheme("x"); A->BinScheme("x")->BuildBins( 6, 0.01, 1, true ); + + // perform the analysis ================================== + A->Execute(); +}; diff --git a/macro/analysis_resolution_SD_Full.C b/macro/analysis_resolution_SD_Full.C new file mode 100644 index 00000000..4823c157 --- /dev/null +++ b/macro/analysis_resolution_SD_Full.C @@ -0,0 +1,41 @@ +R__LOAD_LIBRARY(Largex) + +// make resolution plots +void analysis_resolution_SD_Full( + TString infiles="datarec/full_dis-5x41.config", /* list of input files */ + Double_t eleBeamEn=5, /* electron beam energy [GeV] */ + Double_t ionBeamEn=41, /* ion beam energy [GeV] */ + Double_t crossingAngle=25, /* crossing angle [mrad] */ + TString outfilePrefix="full_dis-5x41" /* output filename prefix*/ +) { + + // setup analysis ======================================== + AnalysisDD4hep *A = new AnalysisDD4hep( + infiles, + eleBeamEn, + ionBeamEn, + crossingAngle, + outfilePrefix + ); + A->NBINS = 5; // use this to set the number of bins along each axis, e.g., z binning (except resolution axes) for each overall bin in e.g. x and Q2 + A->NBINSRES = 100; // use this to set the number of bins along the resolution axes for each overall bin in e.g. x and Q2 + // A->RESHIGH = 1; + // A->RESLOW = -1; + //A->maxEvents = 30000; // use this to limit the number of events + A->SetReconMethod("Ele"); // set reconstruction method + A->AddFinalState("pipTrack"); // pion final state + + // define cuts ==================================== + A->AddBinScheme("w"); A->BinScheme("w")->BuildBin("Min",3.0); // W > 3 GeV + A->AddBinScheme("y"); A->BinScheme("y")->BuildBin("Range",0.01,0.95); // 0.01 < y < 0.95 + A->AddBinScheme("z"); A->BinScheme("z")->BuildBin("Range",0.2,0.9); // 0.2 < z < 0.9 + A->AddBinScheme("xF"); A->BinScheme("xF")->BuildBin("Min",0.0); // xF > 0 + A->AddBinScheme("ptLab"); A->BinScheme("ptLab")->BuildBin("Min",0.1); // pT_lab > 0.1 GeV (tracking limit) + + // set binning scheme ==================================== + A->AddBinScheme("q2"); A->BinScheme("q2")->BuildBins( 4, 1, 100, true ); + A->AddBinScheme("x"); A->BinScheme("x")->BuildBins( 6, 0.01, 1, true ); + + // perform the analysis ================================== + A->Execute(); +}; diff --git a/macro/plotPurities.C b/macro/plotPurities.C new file mode 100644 index 00000000..7bc48df9 --- /dev/null +++ b/macro/plotPurities.C @@ -0,0 +1,586 @@ + + +void plotPurities(){ + // Get purity and efficiency plots from different methods output files + // And plot together + + // TString p1 = "out/1bin/JB_dis-18x275-xm25.canvas.root"; + // TString p2 = "out/1bin/DA_dis-18x275-xm25.canvas.root"; + // TString p3 = "out/1bin/Ele_dis-18x275-xm25.canvas.root"; + + TString p1 = "out/JB_dis-10x275-xm25.canvas.root"; + TString p2 = "out/DA_dis-10x275-xm25.canvas.root"; + TString p3 = "out/Ele_dis-10x275-xm25.canvas.root"; + + TFile *f1 = TFile::Open(p1); + TFile *f2 = TFile::Open(p2); + TFile *f3 = TFile::Open(p3); + + // You should know nbins and limits in x and Q2 already + int nx = 6; + int nq = 4; + double xMin = 1e-2; double xMax = 1; + double qMin = 1; double qMax = 1000; + + TString outName = "dis-10x275-xm25"; + // TString outName = "dis-18x275-xm25"; + const int nNames = 5; + TString histNames[nNames] = {"z_efficiency","z_purity","z_phiH_Res","z_pT_Res","z_z_Res"};//,"z_phiH_Res", + TString labels[nNames] = {"K^{#pm} efficiency","K^{#pm} purity","#phi_{H}","p_{T}","z"};//"#phi_{H}", + TString header = "10x275GeV (0.2 < z < 1.0)"; + // TString header = "18x275GeV (0.2 < z < 1.0)"; + double yMin = -0.1; double yMax=1.0; + + int nbinsz = 1; //Set manually...NOTE TODO + + // borrowed from postprocessor.cxx + + // default values set for nvar1==nvar2 + + // used to be input args + TString var1name="x"; int nvar1=nx; double var1low=xMin; double var1high=xMax; bool var1log=true; + TString var2name="Q2"; int nvar2=nq; double var2low=qMin; double var2high=qMax; bool var2log=true; + bool intlog1=false; bool intlog2=false; bool intgrid1=false; bool intgrid2=false; + + + int canvx = 933;//700; + int canvy = 800;//600;//TODO: check new numbers are better? + double botmargin = 0.2; + double leftmargin = 0.2; + double xaxisy = 0.04; + double xaxisx1 = 0.08; + double xaxisx2 = 0.97; + double yaxisx = 0.04; + double yaxisy1 = 0.085; + double yaxisy2 = 0.97; + if(nvar1 > nvar2){ + // different canvas sizing/axis position for unequal binning + canvx = 1100; + canvy = 700; + xaxisx1 = 0.075; + xaxisx2 = 0.975; + yaxisy1 = 0.08; + } + + TString canvN = "canv_"+outName+"_all__"; + TString histN = "stack_"+outName+"_all__"; + for (int k=0; kSetFrameLineWidth(0); + mainpad->UseCurrentStyle(); + + mainpad->SetFillStyle(4000); + mainpad->Divide(nvar1,nvar2,0,0); + mainpad->Draw(); + TLine * lDIRC = new TLine(6,-1,6,1); + TLine * lDIRClow = new TLine(0.5,-1,0.5,1); + TLine * lmRICH = new TLine(2,-1,2,-4); + TLine * lDRICH = new TLine(2.5,1,2.5,4); + lDIRC->SetLineColor(kRed); + lDIRClow->SetLineColor(kRed); + lmRICH->SetLineColor(kRed); + lDRICH->SetLineColor(kRed); + THStack* histArray[nvar1][nvar2]; + int drawpid = 0; + //outfile->cd("/"); + // canv->Write(); + bool daz = false; + bool dapt = false; + bool daphi = false; + bool elez = false; + bool elept = false; + bool elephi = false; + bool jbz = false; + bool jbpt = false; + bool jbphi = false; + bool purity = false; + bool efficiency = false; + TPad *lgpad = new TPad("lgpad", "lgpad", 0.84, 0.07, 0.98, 0.21); + TLegend *lg = new TLegend(0.01,0.01,0.99,0.99); + lg->SetHeader(header,"C"); + lg->SetTextSize(0.08); + if (nNames>3) lg->SetNColumns(2); + // if (nNames>4) lg->SetNColumns(3); + + for (int i=0; iGet(name)!=nullptr) h1_ = (TH1D*)f1->Get(name); else h1_ = nullptr; + TH1D *h2_; if (f2->Get(name)!=nullptr) h2_ = (TH1D*)f2->Get(name); else h2_ = nullptr; + TH1D *h3_; if (f3->Get(name)!=nullptr) h3_ = (TH1D*)f3->Get(name); else h3_ = nullptr; + + TH1D *h1 = new TH1D("h1","",nbinsz,0,1); + TH1D *h2 = new TH1D("h2","",nbinsz,0,1); + TH1D *h3 = new TH1D("h3","",nbinsz,0,1); + for (int idx=1; idx<=nbinsz; idx++) { + if (h1_!=nullptr && h2_!=nullptr && h3_!=nullptr){ + // h1_->SetBinError(idx,0); h2_->SetBinError(idx,0); h3_->SetBinError(idx,0); + // std::cout<<"\t"<Get(name)<<" "<Get(name)<<" "<Get(name)<GetBinContent(idx)GetBinContent(idx) && h1_->GetBinContent(idx)GetBinContent(idx) && (histNames[k]!="z_purity" && histNames[k]!="z_efficiency")) { + // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetBinContent(idx,h1_->GetBinContent(idx)); + if (histNames[k]=="z_z_Res") h1->SetMarkerStyle(24); + if (histNames[k]=="z_pT_Res") h1->SetMarkerStyle(26); + if (histNames[k]=="z_phiH_Res") h1->SetMarkerStyle(32); + // std::cout<<"\th1: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(8); + h1->SetMarkerSize(1); + h1->GetYaxis()->SetRangeUser(-0.05,1); + hist->Add(h1); + if ( /*((i==4 && j==1) || (i==5 && j==2)) &&*/ ((histNames[k]=="z_z_Res" && !jbz) || (histNames[k]=="z_pT_Res" && !jbpt) || (histNames[k]=="z_phiH_Res" && !jbphi)) ) {//TODO: Find where these actually pop up + + if (histNames[k]=="z_z_Res") jbz=true; + if (histNames[k]=="z_pT_Res") jbpt=true; + if (histNames[k]=="z_phiH_Res") jbphi=true; + std::cout<<"ADDING JB ENTRY "<AddEntry(h1,"JB "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + } + // continue; //IMPORTANT! + } + if (h2_->GetBinContent(idx)GetBinContent(1) && h2_->GetBinContent(idx)GetBinContent(idx) && (histNames[k]!="z_purity" && histNames[k]!="z_efficiency")) { + // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetBinContent(idx,h2_->GetBinContent(idx)); + if (histNames[k]=="z_z_Res") h2->SetMarkerStyle(24); + if (histNames[k]=="z_pT_Res") h2->SetMarkerStyle(26); + if (histNames[k]=="z_phiH_Res") h2->SetMarkerStyle(32); + // std::cout<<"\th2: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(4); + h2->SetMarkerSize(1); + h2->GetYaxis()->SetRangeUser(-0.05,1); + hist->Add(h2); + if ((histNames[k]=="z_z_Res" && !daz) || (histNames[k]=="z_pT_Res" && !dapt) || (histNames[k]=="z_phiH_Res" && !daphi)){//TODO: Find where these actually pop up + if (histNames[k]=="z_z_Res") daz=true; + if (histNames[k]=="z_pT_Res") dapt=true; + if (histNames[k]=="z_phiH_Res") daphi=true; + std::cout<<"ADDING ENTRY "<AddEntry(h2,"DA "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + } + // continue; + } + // else { + // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetBinContent(idx,h3_->GetBinContent(idx)); + if (histNames[k]=="z_z_Res") h3->SetMarkerStyle(24); + if (histNames[k]=="z_pT_Res") h3->SetMarkerStyle(26); + if (histNames[k]=="z_phiH_Res") h3->SetMarkerStyle(32); + if (histNames[k]=="z_purity") h3->SetMarkerStyle(25); + if (histNames[k]=="z_efficiency") h3->SetMarkerStyle(34); + // std::cout<<"\th3: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(2); if (histNames[k]=="z_purity") h3->SetMarkerColor(8); + if (histNames[k]=="z_efficiency") h3->SetMarkerColor(9); + h3->SetMarkerSize(1); + h3->GetYaxis()->SetRangeUser(-0.05,1); + hist->Add(h3); + if ((histNames[k]=="z_z_Res" && !elez) || (histNames[k]=="z_pT_Res" && !elept) || (histNames[k]=="z_phiH_Res" && !elephi) || !purity){//TODO: Find where these actually pop up + if (histNames[k]=="z_z_Res") elez=true; + if (histNames[k]=="z_pT_Res") elept=true; + if (histNames[k]=="z_phiH_Res") elephi=true; + // if (histNames[k]=="z_purity") purity=true; + // if (histNames[k]=="z_efficiency") efficiency=true; + std::cout<<"ADDING ENTRY "<AddEntry(h3,"Ele "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + else if (!purity && histNames[k]=="z_purity") lg->AddEntry(h3,labels[k],"p"); if (histNames[k]=="z_purity") purity=true; + else if (!efficiency && histNames[k]=="z_efficiency") lg->AddEntry(h3,labels[k],"p"); if (histNames[k]=="z_efficiency") efficiency=true; + // } + } + }// if (h1!=nullptr && h2!=nullptr) + + if (h1_==nullptr && h2_!=nullptr && h3_!=nullptr){ + // h2_->SetBinError(idx,0); h3_->SetBinError(idx,0); + // std::cout<<"\t"<Get(name)<<" "<Get(name)<<" "<Get(name)<GetBinContent(1)GetBinContent(1) && h1->GetBinContent(1)GetBinContent(1) && histNames[k]!="z_purity") { + // // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetMarkerStyle(24); + // if (histNames[k]=="z_pT_Res") h1->SetMarkerStyle(26); + // if (histNames[k]=="z_phiH_Res") h1->SetMarkerStyle(32); + // // std::cout<<"\th1: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(8); + // h1->SetMarkerSize(1); + // h1->GetYaxis()->SetRangeUser(-0.05,1); + // hist->Add(h1); + // if ( ((i==4 && j==1) || (i==5 && j==2)) && ((histNames[k]=="z_z_Res" && !jbz) || (histNames[k]=="z_pT_Res" && !jbpt) || (histNames[k]=="z_phiH_Res" && !jbphi)) ) {//TODO: Find where these actually pop up + + // if (histNames[k]=="z_z_Res") jbz=true; + // if (histNames[k]=="z_pT_Res") jbpt=true; + // if (histNames[k]=="z_phiH_Res") jbphi=true; + // std::cout<<"ADDING JB ENTRY "<AddEntry(h1,"JB "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + // } + // // continue; //IMPORTANT! + // } + if (h2_->GetBinContent(idx)GetBinContent(idx) && (histNames[k]!="z_purity" && histNames[k]!="z_efficiency")) { + // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetBinContent(idx,h2_->GetBinContent(idx)); + if (histNames[k]=="z_z_Res") h2->SetMarkerStyle(24); + if (histNames[k]=="z_pT_Res") h2->SetMarkerStyle(26); + if (histNames[k]=="z_phiH_Res") h2->SetMarkerStyle(32); + // std::cout<<"\th2: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(4); + h2->SetMarkerSize(1); + h2->GetYaxis()->SetRangeUser(-0.05,1); + hist->Add(h2); + if ((histNames[k]=="z_z_Res" && !daz) || (histNames[k]=="z_pT_Res" && !dapt) || (histNames[k]=="z_phiH_Res" && !daphi)){//TODO: Find where these actually pop up + if (histNames[k]=="z_z_Res") daz=true; + if (histNames[k]=="z_pT_Res") dapt=true; + if (histNames[k]=="z_phiH_Res") daphi=true; + std::cout<<"ADDING ENTRY "<AddEntry(h2,"DA "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + } + // continue; + } + // else { + // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetBinContent(idx,h3_->GetBinContent(idx)); + if (histNames[k]=="z_z_Res") h3->SetMarkerStyle(24); + if (histNames[k]=="z_pT_Res") h3->SetMarkerStyle(26); + if (histNames[k]=="z_phiH_Res") h3->SetMarkerStyle(32); + if (histNames[k]=="z_purity") h3->SetMarkerStyle(25); + // std::cout<<"\th3: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(2); if (histNames[k]=="z_purity") h3->SetMarkerColor(8); + if (histNames[k]=="z_efficiency") h3->SetMarkerColor(9); + h3->SetMarkerSize(1); + h3->GetYaxis()->SetRangeUser(-0.05,1); + hist->Add(h3); + if ((histNames[k]=="z_z_Res" && !elez) || (histNames[k]=="z_pT_Res" && !elept) || (histNames[k]=="z_phiH_Res" && !elephi) || !purity){//TODO: Find where these actually pop up + if (histNames[k]=="z_z_Res") elez=true; + if (histNames[k]=="z_pT_Res") elept=true; + if (histNames[k]=="z_phiH_Res") elephi=true; + // if (histNames[k]=="z_purity") purity=true; + // if (histNames[k]=="z_efficiency") efficiency=true; + std::cout<<"ADDING ENTRY "<AddEntry(h3,"Ele "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + else if (!purity && histNames[k]=="z_purity") lg->AddEntry(h3,labels[k],"p"); if (histNames[k]=="z_purity") purity=true; + else if (!efficiency && histNames[k]=="z_efficiency") lg->AddEntry(h3,labels[k],"p"); if (histNames[k]=="z_efficiency") efficiency=true; + // } + } + }// if (h2!=nullptr && h3!=nullptr) + + if (h1_!=nullptr && h2_==nullptr && h3_!=nullptr){ + // h1_->SetBinError(1,0); h3_->SetBinError(1,0); + // std::cout<<"\t"<Get(name)<<" "<Get(name)<<" "<Get(name)<GetBinContent(idx)GetBinContent(idx) && (histNames[k]!="z_purity" && histNames[k]!="z_efficiency")) { + h1->SetBinContent(idx,h1_->GetBinContent(idx)); + // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetMarkerStyle(24); + if (histNames[k]=="z_pT_Res") h1->SetMarkerStyle(26); + if (histNames[k]=="z_phiH_Res") h1->SetMarkerStyle(32); + // std::cout<<"\th1: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(8); + h1->SetMarkerSize(1); + h1->GetYaxis()->SetRangeUser(-0.05,1); + hist->Add(h1); + if ( /*((i==4 && j==1) || (i==5 && j==2)) &&*/ ((histNames[k]=="z_z_Res" && !jbz) || (histNames[k]=="z_pT_Res" && !jbpt) || (histNames[k]=="z_phiH_Res" && !jbphi)) ) {//TODO: Find where these actually pop up + + if (histNames[k]=="z_z_Res") jbz=true; + if (histNames[k]=="z_pT_Res") jbpt=true; + if (histNames[k]=="z_phiH_Res") jbphi=true; + std::cout<<"ADDING JB ENTRY "<AddEntry(h1,"JB "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + } + // continue; //IMPORTANT! + } + // if (h2->GetBinContent(1)GetBinContent(1) && histNames[k]!="z_purity") { + // // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetMarkerStyle(24); + // if (histNames[k]=="z_pT_Res") h2->SetMarkerStyle(26); + // if (histNames[k]=="z_phiH_Res") h2->SetMarkerStyle(32); + // // std::cout<<"\th2: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(4); + // h2->SetMarkerSize(1); + // h2->GetYaxis()->SetRangeUser(-0.05,1); + // hist->Add(h2); + // if ((histNames[k]=="z_z_Res" && !daz) || (histNames[k]=="z_pT_Res" && !dapt) || (histNames[k]=="z_phiH_Res" && !daphi)){//TODO: Find where these actually pop up + // if (histNames[k]=="z_z_Res") daz=true; + // if (histNames[k]=="z_pT_Res") dapt=true; + // if (histNames[k]=="z_phiH_Res") daphi=true; + // std::cout<<"ADDING ENTRY "<AddEntry(h2,"DA "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + // } + // // continue; + // } + // else { + // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetBinContent(idx,h3_->GetBinContent(idx)); + if (histNames[k]=="z_z_Res") h3->SetMarkerStyle(24); + if (histNames[k]=="z_pT_Res") h3->SetMarkerStyle(26); + if (histNames[k]=="z_phiH_Res") h3->SetMarkerStyle(32); + if (histNames[k]=="z_purity") h3->SetMarkerStyle(25); + // std::cout<<"\th3: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(2); if (histNames[k]=="z_purity") h3->SetMarkerColor(8); if (histNames[k]=="z_efficiency") h3->SetMarkerColor(8); + h3->SetMarkerSize(1); + h3->GetYaxis()->SetRangeUser(-0.05,1); + hist->Add(h3); + if ((histNames[k]=="z_z_Res" && !elez) || (histNames[k]=="z_pT_Res" && !elept) || (histNames[k]=="z_phiH_Res" && !elephi) || !purity || !efficiency){//TODO: Find where these actually pop up + if (histNames[k]=="z_z_Res") elez=true; + if (histNames[k]=="z_pT_Res") elept=true; + if (histNames[k]=="z_phiH_Res") elephi=true; + + std::cout<<"ADDING ENTRY "<AddEntry(h3,"Ele "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + else if (!purity && histNames[k]=="z_purity") lg->AddEntry(h3,labels[k],"p"); if (histNames[k]=="z_purity") purity=true; + else if (!efficiency && histNames[k]=="z_efficiency") lg->AddEntry(h3,labels[k],"p"); if (histNames[k]=="z_efficiency") efficiency=true; + // } + } + }// if (h1!=nullptr && h3!=nullptr) + } // idx loop + /* + // if (h1!=nullptr && h2!=nullptr && h3==nullptr){ + // // std::cout<<"\t"<Get(name)<<" "<Get(name)<<" "<Get(name)<GetBinContent(1)GetBinContent(1) && histNames[k]!="z_purity") { + // // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetMarkerStyle(24); + // if (histNames[k]=="z_pT_Res") h1->SetMarkerStyle(26); + // if (histNames[k]=="z_phiH_Res") h1->SetMarkerStyle(32); + // // std::cout<<"\th1: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(8); + // h1->SetMarkerSize(1); + // h1->GetYaxis()->SetRangeUser(-0.05,1); + // hist->Add(h1); + // if ( ((i==4 && j==1) || (i==5 && j==2)) && ((histNames[k]=="z_z_Res" && !jbz) || (histNames[k]=="z_pT_Res" && !jbpt) || (histNames[k]=="z_phiH_Res" && !jbphi)) ) {//TODO: Find where these actually pop up + + // if (histNames[k]=="z_z_Res") jbz=true; + // if (histNames[k]=="z_pT_Res") jbpt=true; + // if (histNames[k]=="z_phiH_Res") jbphi=true; + // std::cout<<"ADDING JB ENTRY "<AddEntry(h1,"JB "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + // } + // // continue; //IMPORTANT! + // } + // if (h2->GetBinContent(1)GetBinContent(1) && histNames[k]!="z_purity") { + // // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetMarkerStyle(24); + // if (histNames[k]=="z_pT_Res") h2->SetMarkerStyle(26); + // if (histNames[k]=="z_phiH_Res") h2->SetMarkerStyle(32); + // // std::cout<<"\th2: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(4); + // h2->SetMarkerSize(1); + // h2->GetYaxis()->SetRangeUser(-0.05,1); + // hist->Add(h2); + // if ((histNames[k]=="z_z_Res" && !daz) || (histNames[k]=="z_pT_Res" && !dapt) || (histNames[k]=="z_phiH_Res" && !daphi)){//TODO: Find where these actually pop up + // if (histNames[k]=="z_z_Res") daz=true; + // if (histNames[k]=="z_pT_Res") dapt=true; + // if (histNames[k]=="z_phiH_Res") daphi=true; + // std::cout<<"ADDING ENTRY "<AddEntry(h2,"DA "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + // } + // // continue; + // } + // // // else { + // // // std::cout<<"\tbin content: "<GetBinContent(1)<<" nbins: "<GetNbinsX()<SetMarkerStyle(24); + // // if (histNames[k]=="z_pT_Res") h3->SetMarkerStyle(26); + // // if (histNames[k]=="z_phiH_Res") h3->SetMarkerStyle(32); + // // if (histNames[k]=="z_purity") h3->SetMarkerStyle(25); + // // // std::cout<<"\th3: marker style, color: "<GetMarkerStyle()<<" "<GetMarkerColor()<SetMarkerColor(2); + // // h3->SetMarkerSize(1); + // // h3->GetYaxis()->SetRangeUser(-0.05,1); + // // hist->Add(h3); + // // if ((histNames[k]=="z_z_Res" && !elez) || (histNames[k]=="z_pT_Res" && !elept) || (histNames[k]=="z_phiH_Res" && !elephi) || !purity){//TODO: Find where these actually pop up + // // if (histNames[k]=="z_z_Res") elez=true; + // // if (histNames[k]=="z_pT_Res") elept=true; + // // if (histNames[k]=="z_phiH_Res") elephi=true; + // // if (histNames[k]=="z_purity") purity=true; + // // std::cout<<"ADDING ENTRY "<AddEntry(h3,"Ele "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + // // else lg->AddEntry(h3,labels[k],"p"); + // // // } + // // } + // }// if (h2!=nullptr && h2!=nullptr) + */ + + // if (h1!=nullptr) { + // h1->SetMarkerColor(2); + // h1->SetMarkerSize(1); + // hist->Add(h1); + // if (i==0 && j==0){//TODO: Find where these actually pop up + // lg->AddEntry(h1,"JB "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + // } + // continue; + // } + // if (h2!=nullptr) { + // h2->SetMarkerColor(4); + // h2->SetMarkerSize(1); + // // h1->GetXaxis()->SetNDivision + // hist->Add(h2); + // if (i==1 && j==1){//TODO: Find where these actually pop up + // lg->AddEntry(h2,"DA "+labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + // } + // continue; + // } + } + + mainpad->cd((nvar2-j-1)*nvar1 + i + 1); + gPad->SetLogx(intlog1); + gPad->SetLogy(intlog2); + gPad->SetGridy(intgrid2); + gPad->SetGridx(intgrid1); + // gPad->SetGrid(0,5);//NOTE: ADDED TODO: CHECK + TString drawStr = ""; + switch(1) {//TODO: figure out how to get THStack dimension? //can't use hist->GetHistogram()->GetDimension() + case 1: + if (i==0 || nbinsz!=1) drawStr = "nostack p"; //NOTE: nostackb will just throw an error, don't use. /*"ex0 p nostack"*/ + else drawStr = "nostack p a"; + break; + case 2: + drawStr = "COLZ"; + break; + case 3: + drawStr = "BOX"; + break; + }; + + if( hist->GetNhists() > 0 ) { + hist->Draw(drawStr); + hist->GetHistogram()->GetYaxis()->SetNdivisions(10); + if (nbinsz==1) hist->GetHistogram()->GetXaxis()->SetNdivisions(0); + hist->GetHistogram()->GetYaxis()->SetLabelSize(0.09); + // mainpad->SetGrid(0,0); + // if (i!=0) { + // for (int idx=0; idx<5; idx++){ + // TF1 *f5 = new TF1("f5","0.2",hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + // f5->SetLineStyle(2); + // f5->SetLineColor(1); + // f5->Draw("SAME"); + // TF1 *f2 = new TF1("f2","0.4",hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + // f2->SetLineStyle(2); + // f2->SetLineColor(1); + // f2->Draw("SAME"); + // TF1 *f3 = new TF1("f3","0.6",hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + // f3->SetLineStyle(2); + // f3->SetLineColor(1); + // f3->Draw(); + // TF1 *f4 = new TF1("f4","0.8",hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + // f4->SetLineStyle(2); + // f4->SetLineColor(1); + // f4->Draw("SAME"); + + // } + // } + // hist->GetHistogram()->GetYaxis()->SetRangeUser(0,1); + // hist->Draw(drawStr); + TF1 *f1 = new TF1("f1","0.0",hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + f1->SetLineColor(1); + f1->SetLineStyle(1); + f1->SetLineWidth(1); + f1->Draw("SAME"); + + TF1 *f2 = new TF1("f2","0.2",hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + f2->SetLineColor(1); + f2->SetLineStyle(2); + f2->SetLineWidth(1); + f2->Draw("SAME"); + + TF1 *f3 = new TF1("f3","0.4",hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + f3->SetLineColor(1); + f3->SetLineStyle(2); + f3->SetLineWidth(1); + f3->Draw("SAME"); + + TF1 *f4 = new TF1("f4","0.6",hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + f4->SetLineColor(1); + f4->SetLineStyle(2); + f4->SetLineWidth(1); + f4->Draw("SAME"); + + TF1 *f5 = new TF1("f5","0.8",hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + f5->SetLineColor(1); + f5->SetLineStyle(2); + f5->SetLineWidth(1); + f5->Draw("SAME"); + + TF1 *f6 = new TF1("f6","1.0",hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + f6->SetLineColor(1); + f6->SetLineStyle(2); + f6->SetLineWidth(1); + f6->Draw("SAME"); + + if (i==0 && j==0) { + mainpad->cd(nvar1*nvar2);// Bottom right corner pad + lg->Draw(); + mainpad->cd((nvar2-j-1)*nvar1 + i + 1);// Return to original pad + } + if(drawpid){ + lDIRClow->Draw(); + lDIRC->Draw(); + lmRICH->Draw(); + lDRICH->Draw(); + } + } + + } + } + + // Overall stuff + canv->cd(); + + TPad *newpad1 = new TPad("newpad1","full pad",0,0,1,1); + TPad *newpad2 = new TPad("newpad2","full pad",0,0,1,1); + newpad1->SetFillStyle(4000); + newpad1->Draw(); + newpad2->SetFillStyle(4000); + newpad2->Draw(); + + TString xopt, yopt; + if(var1log) xopt = "GS"; + else xopt = "S"; + if(var2log) yopt = "GS"; + else yopt = "S"; + + TGaxis *xaxis = new TGaxis(xaxisx1,xaxisy,xaxisx2,xaxisy,var1low,var1high,510,xopt); + TGaxis *yaxis = new TGaxis(yaxisx,yaxisy1,yaxisx,yaxisy2,var2low,var2high,510,yopt); + xaxis->SetTitle(var1name); + xaxis->SetName("xaxis"); + xaxis->SetTitleSize(0.02); + xaxis->SetTextFont(40); + xaxis->SetLabelSize(0.02); + xaxis->SetTickSize(0.02); + + yaxis->SetTitle(var2name); + yaxis->SetTitleSize(0.02); + yaxis->SetName("yaxis"); + yaxis->SetTextFont(40); + yaxis->SetLabelSize(0.02); + yaxis->SetTickSize(0.02); + + newpad1->cd(); + yaxis->Draw(); + newpad2->cd(); + xaxis->Draw(); + + // canv->Write(); + canv->Print(canvN+".png"); + canv->Print(canvN+".pdf"); + + // OLD + // // THStack *hs = new THStack("myStack","myStack"); + // // hs->Add(h1); hs->Add(h2); + + // // Make legend + // TString header = "testheader"; + // TLegend *lg = new TLegend(0.80,0.05,0.90,0.15); + // lg->SetHeader(header,"C"); + // lg->SetTextSize(0.15); + // lg->SetNColumns(2); + // lg->AddEntry(h1,"label1","p"); + // lg->AddEntry(h2,"label2","p"); + + // // Draw and save + // TCanvas *c1 = new TCanvas(); + // h1->Draw("BOX"); + // h2->Draw("SAME") + // lg->Draw(); + // c1->Print("myStack.pdf"); + +} \ No newline at end of file diff --git a/macro/postprocess_barak.C b/macro/postprocess_barak.C new file mode 100644 index 00000000..075783b9 --- /dev/null +++ b/macro/postprocess_barak.C @@ -0,0 +1,86 @@ +R__LOAD_LIBRARY(Largex) + +// make resolution plots +// - adapted from `postprocess_pTvsEta.C` +void postprocess_barak( + TString infile="out/barak_dis-18x275-xm25.root" +){ + + gROOT->ProcessLine(".! rm -v out/barak_dis-18x275-xm25.images/*.png"); // cleanup old image files + gROOT->ProcessLine(".! rm -v out/barak_dis-18x275-xm25.images/*.pdf"); // cleanup old image files + + PostProcessor *P = new PostProcessor(infile); + P->Op()->PrintBreadth("HistosDAG Initial Setup"); + + // number of bins in x and Q2 + int nx = P->Op()->GetBinSet("x")->GetNumBins(); + int nq2 = P->Op()->GetBinSet("q2")->GetNumBins(); + + // just counters for filling Histos vector + int xbin = 0; + int q2bin = 0; + + // initialize this 2D vector to be some large size + std::vector> histos_xQ2(30,std::vector(30)); + + auto findxQ2bins = [&histos_xQ2,&P,&xbin,&q2bin,nx,nq2](Histos *H ){ + histos_xQ2[xbin][q2bin] = H; + q2bin++; + if(q2bin == nq2){ + q2bin=0; xbin++; + if(xbin == nx) xbin = 0; + } + }; + + auto drawinxQ2bins = [&histos_xQ2, &P, &nx, &nq2](NodePath *bins){ + TString canvname = "xQ2cov_"; //+ bins->GetVar + for(Node *bin: bins->GetBinNodes()){ + if(bin->GetVarName() == "finalState"){ + canvname+=bin->GetID(); + canvname+="_"; + } + if(bin->GetVarName() == "z"){ + canvname+=bin->GetID(); + canvname+="_"; + } + } + + double xMin = 1e-5; + double xMax = 1; + double q2Min = 1e-1; + double q2Max = 1e4; + + // loop over resolution histograms (see ../src/Analysis.cxx `DefineHist*` calls + // for available histograms, or add your own there) + for( TString histname : {"y_Res"} ) { + P->DrawInBins( + canvname, histos_xQ2, histname, + "x", nx, xMin, xMax, true, + "Q^{2}", nq2, q2Min, q2Max, true + ); + }; + + const int nNames = 1; + double yMin, yMax; yMin = -0.1; yMax = 0.5;//Adjust as needed + TString histNames[nNames] = {"z_y_Res"}; + TString labels[nNames] = {"y"}; + TString header="18x275GeV"; + P->DrawSDInBinsTogether( + canvname, histos_xQ2, header, histNames, labels, nNames, yMin, yMax, + "x", nx, xMin, xMax, true, + "Q^{2}", nq2, q2Min, q2Max, true + ); + }; + + auto beforefunction = [](){ + }; + + P->Op()->Subloop({"x","q2"},beforefunction,drawinxQ2bins); + P->Op()->Payload(findxQ2bins); + + //P->Op()->PrintBreadth("HistosDAG Final Setup"); + + P->Execute(); + + P->Finish(); +}; diff --git a/macro/postprocess_resolution.C b/macro/postprocess_resolution.C new file mode 100644 index 00000000..001d37b4 --- /dev/null +++ b/macro/postprocess_resolution.C @@ -0,0 +1,86 @@ +R__LOAD_LIBRARY(Largex) + +// make resolution plots +// - adapted from `postprocess_pTvsEta.C` +void postprocess_resolution( + TString infile="out/dis-5x41.root" +){ + + gROOT->ProcessLine(".! rm -v out/dis-5x41.images/*.png"); // cleanup old image files + gROOT->ProcessLine(".! rm -v out/dis-5x41.images/*.pdf"); // cleanup old image files + + PostProcessor *P = new PostProcessor(infile); + P->Op()->PrintBreadth("HistosDAG Initial Setup"); + + // number of bins in x and Q2 + int nx = P->Op()->GetBinSet("x")->GetNumBins(); + int nq2 = P->Op()->GetBinSet("q2")->GetNumBins(); + + // just counters for filling Histos vector + int xbin = 0; + int q2bin = 0; + + // initialize this 2D vector to be some large size + std::vector> histos_xQ2(30,std::vector(30)); + + auto findxQ2bins = [&histos_xQ2,&P,&xbin,&q2bin,nx,nq2](Histos *H ){ + histos_xQ2[xbin][q2bin] = H; + q2bin++; + if(q2bin == nq2){ + q2bin=0; xbin++; + if(xbin == nx) xbin = 0; + } + }; + + auto drawinxQ2bins = [&histos_xQ2, &P, &nx, &nq2](NodePath *bins){ + TString canvname = "xQ2cov_"; //+ bins->GetVar + for(Node *bin: bins->GetBinNodes()){ + if(bin->GetVarName() == "finalState"){ + canvname+=bin->GetID(); + canvname+="_"; + } + if(bin->GetVarName() == "z"){ + canvname+=bin->GetID(); + canvname+="_"; + } + } + + double xMin = 1e-4; + double xMax = 1; + double q2Min = 0.99; + double q2Max = 1000; + + // loop over resolution histograms (see ../src/Analysis.cxx `DefineHist*` calls + // for available histograms, or add your own there) + for( TString histname : {"x_Res","y_Res","Q2_Res","phiH_Res","phiS_Res","phiHvsPhiS","z_purity","z_efficiency"} ) { + P->DrawInBins( + canvname, histos_xQ2, histname, + "x", nx, xMin, xMax, true, + "Q^{2}", nq2, q2Min, q2Max, true + ); + }; + + const int nNames = 5; + double yMin, yMax; yMin = -0.1; yMax = 1.0;//Adjust as needed + TString histNames[nNames] = {"z_z_Res","z_pT_Res","z_phiH_Res","z_purity","z_efficiency"}; + TString labels[nNames] = {"z","p_{T}","#phi_{H}","purity","efficiency"}; + TString header="testheader"; + P->DrawSDInBinsTogether( + canvname, histos_xQ2, header, histNames, labels, nNames, yMin, yMax, + "x", nx, xMin, xMax, true, + "Q^{2}", nq2, q2Min, q2Max, true + ); + }; + + auto beforefunction = [](){ + }; + + P->Op()->Subloop({"x","q2"},beforefunction,drawinxQ2bins); + P->Op()->Payload(findxQ2bins); + + //P->Op()->PrintBreadth("HistosDAG Final Setup"); + + P->Execute(); + + P->Finish(); +}; diff --git a/src/Analysis.cxx b/src/Analysis.cxx index fd0fc52c..1c0e8e8d 100644 --- a/src/Analysis.cxx +++ b/src/Analysis.cxx @@ -323,6 +323,44 @@ void Analysis::Prepare() { HS->DefineHist1D("z_Res","z-z^{true}","", NBINS, -1.5, 1.5); HS->DefineHist1D("mX_Res","mX-mX^{true}","GeV", NBINS, -5, 5); HS->DefineHist1D("xF_Res","xF-xF^{true}","", NBINS, -1.5, 1.5); + // OLD + // HS->DefineHist1D("Q2_Res","Q2_{true}-Q2","GeV^{2}", NBINS, -2, 2); + // HS->DefineHist1D("x_Res","x_{true}-x","", NBINS, -2, 2); + // HS->DefineHist1D("y_Res","y_{true}-y","", NBINS, -2, 2); + // HS->DefineHist1D("z_Res","z_{true}-z","", NBINS, -2, 2); + // HS->DefineHist1D("pT_Res","pT_{true}-pT","GeV", NBINS, -2, 2); + // HS->DefineHist1D("phiH_Res","#phi_{h}^{true}-#phi_{h}","", NBINS, -TMath::Pi(), TMath::Pi()); + // HS->DefineHist1D("phiS_Res","#phi_{S}^{true}-#phi_{S}","", NBINS, -TMath::Pi(), TMath::Pi()); + + // resolutions vs. z on x axis. + HS->DefineHist2D("z_Q2_Res","z","","#sigma_{Q2}","", NBINS, 0, 1, NBINSRES, -0.5, 0.5);//TODO: Fill these + HS->DefineHist2D("z_x_Res","z","","#sigma_{x}","", NBINS, 0, 1, NBINSRES, -0.5, 0.5); + HS->DefineHist2D("z_y_Res","z","","#sigma_{y}","", NBINS, 0, 1, NBINSRES, -0.5, 0.5); + HS->DefineHist2D("z_z_Res","z","","#sigma_{z}","", NBINS, 0, 1, NBINSRES, -0.5, 0.5); + HS->DefineHist2D("z_pT_Res","z","","#sigma_{pT}","", NBINS, 0, 1, NBINSRES, -0.5, 0.5); + HS->DefineHist2D("z_phiH_Res","z","","#sigma_{#phi_{h}^{true}}","", NBINS, 0, 1, NBINSRES, -TMath::Pi(), TMath::Pi()); + HS->DefineHist2D("z_phiS_Res","z","","#sigma_{#phi_{S}^{true}}","", NBINS, 0, 1, NBINSRES, -TMath::Pi(), TMath::Pi()); + + // 1D z-binned purity and efficiency + HS->DefineHist1D("z_true","z","", NBINS, 0, 1); + HS->DefineHist1D("z_trueMC","z","", NBINS, 0, 1); + HS->DefineHist1D("z_purity","purity","", NBINS, 0, 1); + HS->DefineHist1D("z_efficiency","efficiency","", NBINS, 0, 1); + + // 1D x-binned purity and efficiency + HS->DefineHist1D("x_true","x","", NBINS, 1e-5, 1, true); + HS->DefineHist1D("x_trueMC","x","", NBINS, 1e-5, 1, true); + HS->DefineHist1D("x_purity","purity","", NBINS, 1e-5, 1, true); + HS->DefineHist1D("x_efficiency","efficiency","", NBINS, 1e-5, 1, true); + + // // 2D Q2 vs. x binned resolutions + // HS->DefineHist1D("x_Res","x-x_{true}","", NBINS, -0.5, 0.5); + // HS->DefineHist1D("y_Res","y-y_{true}","", NBINS, -0.2, 0.2); + // HS->DefineHist1D("Q2_Res","Q2-Q2_{true}","GeV", NBINS, -0.5, 0.5); + // HS->DefineHist1D("phiH_Res","#phi_{h}-#phi_{h}^{true}","", NBINS, -TMath::Pi(), TMath::Pi()); + // HS->DefineHist1D("phiS_Res","#phi_{S}-#phi_{S}^{true}","", NBINS, -0.1*TMath::Pi(), 0.1*TMath::Pi()); + // HS->DefineHist1D("pT_Res","pT-pT^{true}","GeV", NBINS, -1.5, 1.5); + HS->DefineHist2D("Q2vsXtrue","x","Q^{2}","","GeV^{2}", 20,1e-4,1, 10,1,1e4, @@ -365,7 +403,6 @@ void Analysis::Prepare() { }); HD->ExecuteAndClearOps(); - // initialize total weights wTrackTotal = 0.; wJetTotal = 0.; @@ -425,7 +462,6 @@ Int_t Analysis::GetEventQ2Idx(Double_t Q2, Int_t guess) { } } - // finish the analysis //----------------------------------- void Analysis::Finish() { @@ -442,18 +478,25 @@ void Analysis::Finish() { // calculate cross sections, and print yields HD->Initial([this](){ cout << sep << endl << "Histogram Entries:" << endl; }); HD->Final([this](){ cout << sep << endl; }); - HD->Payload([&lumi](Histos *H){ + HD->Payload([&lumi,this](Histos *H){ cout << H->GetSetTitle() << " ::: " << H->Hist("Q2vsX")->GetEntries() << endl; // calculate cross sections H->Hist("Q_xsec")->Scale(1./lumi); // TODO: generalize (`if (name contains "xsec") ...`) - // divide resolution plots by true counts per x-Q2 bin - H->Hist("Q2vsXpurity")->Divide(H->Hist("Q2vsXtrue")); - H->Hist("Q2vsX_zres")->Divide(H->Hist("Q2vsXtrue")); - H->Hist("Q2vsX_pTres")->Divide(H->Hist("Q2vsXtrue")); - H->Hist("Q2vsX_phiHres")->Divide(H->Hist("Q2vsXtrue")); + + + // Convert to contamination plot since the y scale is better for plotting with stddevs + // H->Hist("z_purity")->Add(H->Hist("z_true"),-1); + H->Hist("z_purity")->Divide(H->Hist("z_true")); + H->Hist("z_efficiency")->Divide(H->Hist("z_trueMC")); + H->Hist("x_purity")->Divide(H->Hist("x_true")); + H->Hist("x_efficiency")->Divide(H->Hist("x_trueMC")); + // TF1 *f1 = new TF1("f1","1",0,1); //TODO: Set limits automatically + // H->Hist("z_purity")->Multiply(f1,-1); + }); + HD->ExecuteAndClearOps(); // write histograms @@ -653,6 +696,19 @@ void Analysis::FillHistosTracks() { H->Hist("Nu_Res")->Fill( kin->Nu - kinTrue->Nu, wTrack ); H->Hist("phiH_Res")->Fill( Kinematics::AdjAngle(kin->phiH - kinTrue->phiH), wTrack ); H->Hist("phiS_Res")->Fill( Kinematics::AdjAngle(kin->phiS - kinTrue->phiS), wTrack ); + + // z binned resolutions + if(kinTrue->Q2!=0) dynamic_cast(H->Hist("z_Q2_Res"))->Fill( kinTrue->z, (kinTrue->Q2 - kin->Q2)/kinTrue->Q2, wTrack ); + if(kinTrue->z!=0) dynamic_cast(H->Hist("z_x_Res"))->Fill( kinTrue->z, (kinTrue->x - kin->x)/kinTrue->x, wTrack ); + if(kinTrue->y!=0) dynamic_cast(H->Hist("z_y_Res"))->Fill( kinTrue->z, (kinTrue->y - kin->y)/kinTrue->y, wTrack ); + if(kinTrue->z!=0) dynamic_cast(H->Hist("z_z_Res"))->Fill( kinTrue->z, (kinTrue->z - kin->z)/kinTrue->z, wTrack ); + if(kinTrue->pT!=0) dynamic_cast(H->Hist("z_pT_Res"))->Fill( kinTrue->z, (kinTrue->pT - kin->pT)/kinTrue->pT, wTrack ); + dynamic_cast(H->Hist("z_phiH_Res"))->Fill( kinTrue->z, Kinematics::AdjAngle(kin->phiH - kinTrue->phiH), wTrack ); + dynamic_cast(H->Hist("z_phiS_Res"))->Fill( kinTrue->z, Kinematics::AdjAngle(kin->phiS - kinTrue->phiS), wTrack ); + + // purities + // H->Hist("z_true")->Fill(kinTrue->z, wTrack ); + // if( (H->Hist("z_true"))->FindBin(kinTrue->z) == (H->Hist("z_true"))->FindBin(kin->z) ) H->Hist("z_purity")->Fill(kin->z,wTrack); H->Hist("pT_Res")->Fill( kin->pT - kinTrue->pT, wTrack ); H->Hist("z_Res")->Fill( kin->z - kinTrue->z, wTrack ); H->Hist("mX_Res")->Fill( kin->mX - kinTrue->mX, wTrack ); @@ -678,6 +734,70 @@ void Analysis::FillHistosTracks() { HD->ExecuteOps(true); }; +// purity +void Analysis::FillHistosPurity(bool recMatch, bool mcMatch) { + + // add kinematic values to `valueMap` + valueMap.clear(); + activeEvent = false; + /* DIS */ + valueMap.insert(std::pair( "x", kin->x )); + valueMap.insert(std::pair( "q2", kin->Q2 )); + valueMap.insert(std::pair( "y", kin->y )); + valueMap.insert(std::pair( "z", kin->z )); + + // check bins + // - activates HistosDAG bin nodes which contain this track + // - sets `activeEvent` if there is at least one multidimensional bin to fill + HD->TraverseBreadth(CheckBin()); + if(!activeEvent) return; + + // fill histograms, for activated bins only + HD->Payload([this,recMatch,mcMatch](Histos *H){ + if (recMatch) H->Hist("z_true")->Fill(kinTrue->z, wTrack ); + if (mcMatch) H->Hist("z_purity")->Fill(kinTrue->z,wTrack); + if (recMatch) H->Hist("x_true")->Fill(kinTrue->x, wTrack ); + if (mcMatch) H->Hist("x_purity")->Fill(kinTrue->x,wTrack); + }); + // execute the payload + // - save time and don't call `ClearOps` (next loop will overwrite lambda) + // - called with `activeNodesOnly==true` since we only want to fill bins associated + // with this jet + HD->ExecuteOps(true); +}; + +// purity +void Analysis::FillHistosEfficiency(bool noMatch, bool mcMatch) { + + // add kinematic values to `valueMap` + valueMap.clear(); + activeEvent = false; + /* DIS */ + valueMap.insert(std::pair( "x", kin->x )); + valueMap.insert(std::pair( "q2", kin->Q2 )); + valueMap.insert(std::pair( "y", kin->y )); + valueMap.insert(std::pair( "z", kin->z )); + + // check bins + // - activates HistosDAG bin nodes which contain this track + // - sets `activeEvent` if there is at least one multidimensional bin to fill + HD->TraverseBreadth(CheckBin()); + if(!activeEvent) return; + + // fill histograms, for activated bins only + HD->Payload([this,noMatch,mcMatch](Histos *H){ + if (noMatch) H->Hist("z_trueMC")->Fill(kinTrue->z, wTrack ); + if (mcMatch) H->Hist("z_efficiency")->Fill(kinTrue->z,wTrack); + if (noMatch) H->Hist("x_trueMC")->Fill(kinTrue->x, wTrack ); + if (mcMatch) H->Hist("x_efficiency")->Fill(kinTrue->x,wTrack); + }); + // execute the payload + // - save time and don't call `ClearOps` (next loop will overwrite lambda) + // - called with `activeNodesOnly==true` since we only want to fill bins associated + // with this jet + HD->ExecuteOps(true); +}; + // jets void Analysis::FillHistosJets() { @@ -721,9 +841,6 @@ void Analysis::FillHistosJets() { HD->ExecuteOps(true); }; - - - // destructor Analysis::~Analysis() { }; diff --git a/src/Analysis.h b/src/Analysis.h index f168491f..3b0e84cd 100644 --- a/src/Analysis.h +++ b/src/Analysis.h @@ -17,6 +17,7 @@ #include "TClonesArray.h" #include "TFile.h" #include "TRegexp.h" +#include "TF1.h" // largex-eic #include "Histos.h" @@ -46,7 +47,8 @@ class Analysis : public TNamed ~Analysis(); // number of bins for histograms - const Int_t NBINS = 50; + Int_t NBINS = 50; + Int_t NBINSRES = 50; const Int_t NBINS_FULL = 10; // bin schemes @@ -98,6 +100,8 @@ class Analysis : public TNamed // FillHistos methods: fill histograms void FillHistosTracks(); + void FillHistosPurity(bool recMatch, bool mcMatch); + void FillHistosEfficiency(bool noMatch, bool mcMatch); void FillHistosJets(); // lambda to check which bins an observable is in, during DAG breadth diff --git a/src/AnalysisDD4hep.cxx b/src/AnalysisDD4hep.cxx index 30e4a26c..a1362901 100644 --- a/src/AnalysisDD4hep.cxx +++ b/src/AnalysisDD4hep.cxx @@ -266,6 +266,8 @@ void AnalysisDD4hep::process_event() { if(mcid_ == imc.mcID) { + bool mcMatch=(pid_==imc.pid); + FillHistosPurity(false,mcMatch); kinTrue->vecHadron = imc.vecPart; break; } @@ -283,6 +285,7 @@ void AnalysisDD4hep::process_event() if( deta < mineta ) { mineta = deta; + FillHistosPurity(pid_,mcpart[imc].pid); kinTrue->vecHadron = mcpart[imc].vecPart; } } diff --git a/src/AnalysisDelphes.cxx b/src/AnalysisDelphes.cxx index e10a651f..de28d3b4 100644 --- a/src/AnalysisDelphes.cxx +++ b/src/AnalysisDelphes.cxx @@ -26,7 +26,6 @@ AnalysisDelphes::AnalysisDelphes( /* ... none defined yet ... */ }; - //============================================= // perform the analysis //============================================= @@ -49,7 +48,7 @@ void AnalysisDelphes::Execute() { // calculate cross section if(maxEvents>0) ENT = maxEvents; // limiter - // branch iterators + // branch iterators (NEW) TObjArrayIter itTrack(tr->UseBranch("Track")); TObjArrayIter itElectron(tr->UseBranch("Electron")); TObjArrayIter itParticle(tr->UseBranch("Particle")); @@ -164,7 +163,7 @@ void AnalysisDelphes::Execute() { // get vector of jets // TODO: should this have an option for clustering method? kin->GetJets(itEFlowTrack, itEFlowPhoton, itEFlowNeutralHadron, itParticle); - + // track loop - - - - - - - - - - - - - - - - - - - - - - - - - - - - itTrack.Reset(); while(Track *trk = (Track*) itTrack()) { @@ -173,7 +172,12 @@ void AnalysisDelphes::Execute() { // final state cut // - check PID, to see if it's a final state we're interested in for // histograms; if not, proceed to next track - // pid = trk->PID; //NOTE: trk->PID is currently not smeared so it just returns the truth-level PID + + // MC Truth PID + int mcpid = trk->PID; //NOTE: trk->PID is currently not smeared so it just returns the truth-level PID + auto kvMC = PIDtoFinalState.find(mcpid); + + // Reconstructed PID pid = kin->getTrackPID( // get smeared PID trk, itpfRICHTrack, @@ -182,11 +186,16 @@ void AnalysisDelphes::Execute() { itdualRICHagTrack, itdualRICHcfTrack ); auto kv = PIDtoFinalState.find(pid); - if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue; - if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue; - // get parent particle, to check if pion is from vector meson + // if(kv!=PIDtoFinalState.end()) { finalStateID = kv->second; + // if(activeFinalStates.find(finalStateID)!=activeFinalStates.end()) { + if (true) { finalStateID = "pipTrack"; if (pid==321 || pid==-321 || pid==211 || pid==-211) { //DEBUGGING: This is just a temporary fix for not being able to select multiple final states in the same bin. + + // Get total # of final state particles correctly identified in reconstruction GenParticle *trkParticle = (GenParticle*)trk->Particle.GetObject(); + + // get parent particle, to check if pion is from vector meson + // GenParticle *trkParticle = (GenParticle*)trk->Particle.GetObject(); TObjArray *brParticle = (TObjArray*)itParticle.GetCollection(); GenParticle *parentParticle = (GenParticle*)brParticle->At(trkParticle->M1); int parentPID = (parentParticle->PID); // TODO: this is not used yet... @@ -221,16 +230,71 @@ void AnalysisDelphes::Execute() { wTrack = Q2weightFactor * weight->GetWeight(*kinTrue); wTrackTotal += wTrack; + // Get total # of final state particles identified in selected final state + if (pid==321 || pid==-321) FillHistosPurity(true,false); + + // And number correctly identified + mcpid = trkParticle->PID; + if ((pid==321 || pid==-321) && pid==mcpid) { + FillHistosPurity(false,true); + FillHistosEfficiency(true,true); + } + // fill track histograms in activated bins - FillHistosTracks(); + if (pid==211 || pid==-211) FillHistosTracks(); + + // fill simple tree + // - not binned + // - `activeEvent` is only true if at least one bin gets filled for this track + if( writeSimpleTree && activeEvent ) ST->FillTree(wTrack); + + // tests + //kin->ValidateHeadOnFrame(); + + } // if(kv!=PIDtoFinalState.end()) + } // if(activeFinalStates.find(finalStateID)!=activeFinalStates.end()) + if (pid!=mcpid) { + // if(kvMC!=PIDtoFinalState.end()) { finalStateID = kv->second; + // if(activeFinalStates.find(finalStateID)!=activeFinalStates.end()) { + if (true) { finalStateID = "pipTrack"; if (mcpid==321 || mcpid==-321 || mcpid==211 || mcpid==-211) { + + // calculate hadron kinematics + GenParticle* trkPart = (GenParticle*)trk->Particle.GetObject(); + kinTrue->hadPID = mcpid; + kinTrue->vecHadron.SetPtEtaPhiM( + trkPart->PT, + trkPart->Eta, + trkPart->Phi, + trkPart->Mass /* TODO: do we use track mass here ?? */ + ); + + kinTrue->CalculateHadronKinematics(); + + // asymmetry injection + //kin->InjectFakeAsymmetry(); // sets tSpin, based on reconstructed kinematics + //kinTrue->InjectFakeAsymmetry(); // sets tSpin, based on generated kinematics + //kin->tSpin = kinTrue->tSpin; // copy to "reconstructed" tSpin + + // Get index of file that the event comes from. + Double_t Q2weightFactor = GetEventQ2Weight(kinTrue->Q2, inLookup[chain->GetTreeNumber()]); + wTrack = Q2weightFactor * weight->GetWeight(*kinTrue); + wTrackTotal += wTrack; + + // Get total # of final state particles identified in selected final state + if (mcpid==321 || mcpid==-321) FillHistosEfficiency(true,false); //NOTE: DEBUGGING For now just quick fix since can't combine multiple final states in a bin. // fill simple tree // - not binned // - `activeEvent` is only true if at least one bin gets filled for this track + // - TODO [critical]: add a `finalState` cut (also needed in AnalysisDD4hep) if( writeSimpleTree && activeEvent ) ST->FillTree(wTrack); // tests //kin->ValidateHeadOnFrame(); + + } // if (activeFinalStates.find(finalStateID)!=activeFinalStates.end()) + } // if (kv!=PIDtoFinalState.end()) + } // if (pid!=mcpid) }; // end track loop @@ -273,7 +337,6 @@ void AnalysisDelphes::Execute() { cout << "end event loop" << endl; // event loop end ========================================================= - // finish execution Finish(); cout << "DEBUG PID: nSmeared=" << kin->countPIDsmeared << " nNotSmeared=" << kin->countPIDtrue << endl; diff --git a/src/PostProcessor.cxx b/src/PostProcessor.cxx index c3c62074..b2db7d0c 100644 --- a/src/PostProcessor.cxx +++ b/src/PostProcessor.cxx @@ -236,6 +236,46 @@ void PostProcessor::DrawSingle(Histos *H, TString histName, TString drawFormat) canv->Print(pngDir+"/"+canvN+".png"); }; +/* ALGORITHM: draw a summary histogram to a canvas, and write it + * - `histName` is the name of the histogram given in Histos + * - `drawFormat` is the formatting string passed to TH1::Draw() + */ +void PostProcessor::DrawSummary(Histos *H, TString histName, TString drawFormat) { + cout << "draw summary plot " << histName << "..." << endl; + TH1 *hist = H->Hist(histName); + if(hist==nullptr) { + cerr << "ERROR: cannot find histogram " << histName << endl; + return; + }; + + // Repurposed from deprecated code in OLD VERSION of DrawSingle() below. + TH1 *histClone = (TH1*) hist->Clone(); + TString canvN = "summaryCanv_"+histName; + histClone->SetLineColor(nsumSetMarkerColor(nsumSetMarkerStyle(nsumSetGrid(1,1); + summaryCanv->SetLogx(H->GetHistConfig(histName)->logx); + summaryCanv->SetLogy(H->GetHistConfig(histName)->logy); + summaryCanv->SetLogz(H->GetHistConfig(histName)->logz); + summaryCanv->SetBottomMargin(0.15); + summaryCanv->SetLeftMargin(0.15); + }; + summaryCanv->cd(); + // std::cout << "histClone max = " << histClone->GetBinContent(histClone->GetMaximumBin()) << std::endl; //DEBUGGING + histClone->Draw(drawFormat+(nsum>0?" SAME":"")); + summaryCanv->Print(pngDir+"/"+canvN+".png"); + nsum++; + +}; + // OLD VERSION: void PostProcessor::DrawSingle(TString histSet, TString histName) { @@ -373,6 +413,8 @@ void PostProcessor::DrawInBins( yaxisy1 = 0.08; }; + double yMin = 0; double yMax = 1;//TODO CHECK THIS WORKS + TString canvN = "canv_"+outName+"_"+histName; TCanvas *canv = new TCanvas(canvN,canvN, canvx, canvy); TPad *mainpad = new TPad("mainpad", "mainpad", 0.07, 0.07, 0.98, 0.98); @@ -397,13 +439,29 @@ void PostProcessor::DrawInBins( for(int j = 0; j < nvar2; j++){ //Histos *H = (Histos*) infile->Get(histList[i][j]); Histos *H = histList[i][j]; + // INTRODUCE LOOP OVER HISTNAMES HERE -> TURN INTO MULTIGRAPHS AND SHOW PERSPECTIVE WISE TH1 *hist = H->Hist(histName); histArray[i][j] = hist; hist->SetTitle(""); - //hist->GetXaxis()->SetTitle(""); - //hist->GetYaxis()->SetTitle(""); - //hist->GetXaxis()->SetLabelSize(0); - //hist->GetYaxis()->SetLabelSize(0); + // //hist->GetXaxis()->SetTitle(""); + // //hist->GetYaxis()->SetTitle(""); + // //hist->GetXaxis()->SetLabelSize(0); + // //hist->GetYaxis()->SetLabelSize(0); + // hist->GetXaxis()->SetTitleSize(0.1); + // hist->GetXaxis()->SetTitleOffset(0.5); + // hist->GetXaxis()->SetNdivisions(8); + // hist->GetXaxis()->SetLabelSize(0.06); + // hist->GetXaxis()->CenterTitle(); + // hist->GetXaxis()->SetLabelOffset(0.02); + // hist->GetYaxis()->SetRangeUser(0,0.05);//TODO: CHECK THIS IS REASONABLE ALSO WHAT ABOUT ERROR BARS??? + // hist->GetYaxis()->SetNdivisions(8); + // hist->GetYaxis()->SetLabelSize(0.06); + // hist->GetYaxis()->SetLabelOffset(0.02); + // // for(int k = 0; k < i; k++){ + // // for(int l = 0; l < j; k++){ + // // histArray[k][l]->GetYaxis()->SetRangeUser(TMath::Min(yMin,hist->GetYaxis()->GetMinimum(),TMath::Max(yMax,hist->GetYaxis()->GetMaximum()))); + // // } + // // } mainpad->cd((nvar2-j-1)*nvar1 + i + 1); gPad->SetLogx(H->GetHistConfig(histName)->logx); @@ -424,8 +482,10 @@ void PostProcessor::DrawInBins( break; }; //hist->Write(); - if( hist->GetEntries() > 0 ) { + if( hist->GetEntries() > 0 ) { //TODO come back and do this after setting up all histograms??? hist->Draw(drawStr); + // TF1 *zeroF = new TF1("zeroF","0",0,1); //TODO: Add optional switch for this? + // zeroF->Draw("SAME"); //TODO: Added zero line if(drawpid){ lDIRClow->Draw(); lDIRC->Draw(); @@ -494,6 +554,265 @@ void PostProcessor::DrawInBins( }; //========================================================================= +/* Convert 2D histogram to 1D histogram of stddevs from y distribution. + * Fits slices along x axis with Gaussian. + */ + +TH1D *PostProcessor::GetSDs(TH2D* fitHist){ + + int nbins = fitHist->GetNbinsX(); + double high = fitHist->GetXaxis()->GetXmax(); + double low = fitHist->GetXaxis()->GetXmin(); + double highY = fitHist->GetYaxis()->GetXmax(); + double lowY = fitHist->GetYaxis()->GetXmin(); + TString name = fitHist->GetName(); + TString title = fitHist->GetTitle(); + TString xtitle = fitHist->GetXaxis()->GetTitle(); + TH1D* subHist = new TH1D(name,title,nbins,low,high); + subHist->GetXaxis()->SetTitle(xtitle); + + // Loop x-axis bins and fit y profiles to get stdev and errors + for (int bin=1; bin<=nbins; bin++) { //NOTE: Bin #'s begin at 1! + + // Define fit function //NOTE: Keep this definition, do not change syntax or use lambda! It will not work, still not sure why... + TF1 *func = new TF1("func","[0]*(1/([1]*TMath::Sqrt(2*TMath::Pi())))*TMath::Exp(-0.5*(x-[2])*(x-[2])/[1]/[1])",lowY,highY); + TH1D *h = new TH1D("h","h",fitHist->GetNbinsY(),lowY,highY); + for (int i=1; i<=h->GetNbinsX(); i++) { //NOTE: Bin #'s start at 1! + h->SetBinContent(i,fitHist->GetBinContent(bin,i)); + } + if (h->GetEntries() < 10 || h->GetMaximum()<1 ) continue; + func->SetParameters(h->GetMaximum()*h->GetStdDev(),h->GetStdDev(),h->GetMean()); + TFitResultPtr fr = h->Fit("func","NS","",lowY,highY); //IMPORTANT: N option keeps fit results from being plotted, otherwise you change the current canvas. + if (fr->IsEmpty() || !(fr->IsValid())) { continue; } + // TMatrixDSym *covMat = new TMatrixDSym(fr->GetCovarianceMatrix()); + + // Gaussian fit parameters + double N0 = func->GetParameter(0); + double sigma = func->GetParameter(1); + double mu = func->GetParameter(2); + + // Gaussian fit errors + double EN0 = func->GetParError(0); + double Esigma = func->GetParError(1); + double Emu = func->GetParError(2); + double chi2 = func->GetChisquare(); + double ndf = func->GetNDF(); + + // //DEBUGGING + // std::cout<<"\tINFO bin["<SetFillStyle(4000); + mainpad->Divide(nvar1,nvar2,0,0); + mainpad->Draw(); + TLine * lDIRC = new TLine(6,-1,6,1); + TLine * lDIRClow = new TLine(0.5,-1,0.5,1); + TLine * lmRICH = new TLine(2,-1,2,-4); + TLine * lDRICH = new TLine(2.5,1,2.5,4); + lDIRC->SetLineColor(kRed); + lDIRClow->SetLineColor(kRed); + lmRICH->SetLineColor(kRed); + lDRICH->SetLineColor(kRed); + THStack* histArray[nvar1][nvar2]; + int drawpid = 0; + outfile->cd("/"); + canv->Write(); + + // get histograms from Histos 2D vector + for(int i = 0; i < nvar1; i++){ + for(int j = 0; j < nvar2; j++){ + //Histos *H = (Histos*) infile->Get(histList[i][j]); + Histos *H = histList[i][j]; + THStack *hist = new THStack(); + TLegend *lg = new TLegend(0.05,0.05,0.95,0.95); + lg->SetHeader(header,"C"); + lg->SetTextSize(0.15); + if (nNames>3) lg->SetNColumns(2); + + for (int k=0; kGetXaxis()->SetTitle(""); + //subHist->GetYaxis()->SetTitle(""); + //subHist->GetXaxis()->SetLabelSize(0); + //subHist->GetYaxis()->SetLabelSize(0); + TH1D *subHist; + if (histNames[k]=="z_purity" || histNames[k]=="z_efficiency" || histNames[k]=="x_purity" || histNames[k]=="x_efficiency") subHist = (TH1D*)H->Hist(histNames[k])->Clone(); + else { + TH2D *fitHist = (TH2D*)H->Hist(histNames[k])->Clone(); + if ( fitHist->GetEntries() < 10 ) continue; //NOTE: Filter out low filled hists that can't get good fits. + fitHist->SetTitle(""); + //DEBUGGING + TString fithistname; fithistname.Form("fithist_"+histNames[k]+"_bin_%d_%d",i,j); + fitHist->SetName(fithistname); fitHist->Write(); + //DEBUGGING + subHist = this->GetSDs(fitHist); + } + + subHist->GetXaxis()->SetTitleSize(0.1); + subHist->GetXaxis()->SetTitleOffset(0.5); + subHist->GetXaxis()->SetNdivisions(8); + subHist->GetXaxis()->SetLabelSize(0.1); + subHist->GetXaxis()->CenterTitle(); + subHist->GetXaxis()->SetLabelOffset(0.02); + subHist->GetYaxis()->SetRangeUser(yMin,yMax);//TODO: CHECK THIS IS REASONABLE ALSO WHAT ABOUT ERROR BARS?? + subHist->GetYaxis()->SetNdivisions(8); + subHist->GetYaxis()->SetLabelSize(0.1); + subHist->GetYaxis()->SetLabelOffset(0.02); + + subHist->SetTitle(histNames[k]); + subHist->SetMarkerStyle(k==0 ? 32 : k+25); + subHist->SetMarkerColor(k+2); + if (k+2>=5) subHist->SetMarkerColor(k+3); //NOTE: 5 is yellow: very hard to see. + subHist->SetMarkerSize(0.5);//NOTE: Remember these will be small plots so keep the binning small and the markers big + if ( subHist->GetEntries()>0 || ((histNames[k]=="z_purity" || histNames[k]=="z_efficiency" || histNames[k]=="x_purity" || histNames[k]=="x_efficiency") && H->Hist("z_z_Res")->GetMaximum()!=0)) { + hist->Add(subHist); + TString myname; myname.Form("hist__"+histNames[k]+"__bin_%d_%d",i,j); + subHist->SetName(myname); + subHist->Write(); + + if (i==0 && j==0){ + lg->AddEntry(subHist,labels[k],"p");//NOTE: Only grabs hists that are in 0,0 bin + } + + } + } + histArray[i][j] = hist; + + mainpad->cd((nvar2-j-1)*nvar1 + i + 1); + gPad->SetLogx(intlog1); + gPad->SetLogy(intlog2); + gPad->SetGridy(intgrid2); + gPad->SetGridx(intgrid1); + TString drawStr = ""; + switch(1) {//TODO: figure out how to get THStack dimension? //can't use hist->GetHistogram()->GetDimension() + case 1: + drawStr = "hist p nostack"; //NOTE: nostackb will just throw an error, don't use. /*"ex0 p nostack"*/ + break; + case 2: + drawStr = "COLZ"; + break; + case 3: + drawStr = "BOX"; + break; + }; + hist->Write(); + if( hist->GetNhists() > 0 ) { + hist->Draw(drawStr); + TF1 *f1 = new TF1("f1","0",hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + f1->SetLineColor(1); + f1->SetLineWidth(1); + f1->Draw("SAME"); + if (i==0 && j==0) { + mainpad->cd(nvar1*nvar2);// Bottom right corner pad + lg->Draw(); + mainpad->cd((nvar2-j-1)*nvar1 + i + 1);// Return to original pad + } + if(drawpid){ + lDIRClow->Draw(); + lDIRC->Draw(); + lmRICH->Draw(); + lDRICH->Draw(); + } + } + }; + }; + canv->cd(); + + TPad *newpad1 = new TPad("newpad1","full pad",0,0,1,1); + TPad *newpad2 = new TPad("newpad2","full pad",0,0,1,1); + newpad1->SetFillStyle(4000); + newpad1->Draw(); + newpad2->SetFillStyle(4000); + newpad2->Draw(); + + TString xopt, yopt; + if(var1log) xopt = "GS"; + else xopt = "S"; + if(var2log) yopt = "GS"; + else yopt = "S"; + + TGaxis *xaxis = new TGaxis(xaxisx1,xaxisy,xaxisx2,xaxisy,var1low,var1high,510,xopt); + TGaxis *yaxis = new TGaxis(yaxisx,yaxisy1,yaxisx,yaxisy2,var2low,var2high,510,yopt); + xaxis->SetTitle(var1name); + xaxis->SetName("xaxis"); + xaxis->SetTitleSize(0.02); + xaxis->SetTextFont(40); + xaxis->SetLabelSize(0.02); + xaxis->SetTickSize(0.02); + + yaxis->SetTitle(var2name); + yaxis->SetTitleSize(0.02); + yaxis->SetName("yaxis"); + yaxis->SetTextFont(40); + yaxis->SetLabelSize(0.02); + yaxis->SetTickSize(0.02); + + newpad1->cd(); + yaxis->Draw(); + newpad2->cd(); + xaxis->Draw(); + + // canv->Write(); + canv->Print(pngDir+"/"+canvN+".png"); + canv->Print(pngDir+"/"+canvN+".pdf"); + outfile->cd("/"); + canv->Write(); + for(int i = 0; i Write(); + } + } +}; + +// ========================================================================= /* ALGORITHM: draw a ratio of all 1D histograms in the specified histogram set * - the ratio will be of `numerSet` over `denomSet` diff --git a/src/PostProcessor.h b/src/PostProcessor.h index b58b2b6c..5e145040 100644 --- a/src/PostProcessor.h +++ b/src/PostProcessor.h @@ -3,6 +3,9 @@ #include #include +#include +#include +#include // root #include "TFile.h" @@ -12,6 +15,12 @@ #include "TROOT.h" #include "TStyle.h" #include "TGaxis.h" +#include "THStack.h" +#include "TLegend.h" +#include "TF1.h" +#include "TFitResult.h" +#include "TFitResultPtr.h" +#include "TMatrixDSym.h" // largex-eic #include "Histos.h" @@ -34,7 +43,6 @@ class PostProcessor : public TNamed const Int_t dimy=700; static const int nsumMax=3; // number of summary plots with formatting - // DAG interfaces: HistosDAG *GetHistosDAG() { return HD; }; HistosDAG *Op() { return GetHistosDAG(); }; // syntactic sugar @@ -62,6 +70,7 @@ class PostProcessor : public TNamed void DumpHist(TString datFile, TString histSet, TString varName); void DumpAve(TString datFile, Histos *H, TString cutName); void DrawSingle(Histos *H, TString histName, TString drawFormat=""); + void DrawSummary(Histos *H, TString histName, TString drawFormat=""); void DrawSingle(TString histSet, TString histName); void DrawRatios( TString outName, Histos *numerSet, Histos *denomSet, Bool_t plotRatioOnly=false @@ -74,6 +83,16 @@ class PostProcessor : public TNamed bool intgrid1=false, bool intgrid2=false ); + TH1D *GetSDs(TH2D* fitHist); + + void DrawSDInBinsTogether( + TString outName, + std::vector>& histList, TString header, TString histNames[], TString labels[], int nNames, double yMin, double yMax, + TString var1name, int nvar1, double var1low, double var1high, bool var1log, + TString var2name, int nvar2, double var2low, double var2high, bool var2log, + bool intlog1=false, bool intlog2=false, bool intgrid1=false, bool intgrid2=false + ); + // algorithm finish methods; to be called after loops void FinishDumpAve(TString datFile); void FinishDrawRatios(TString summaryDir); diff --git a/submit.sh b/submit.sh new file mode 100755 index 00000000..78cecc0e --- /dev/null +++ b/submit.sh @@ -0,0 +1,13 @@ +#!/bin/bash + +#SBATCH --job-name=athenaSIDIS +#SBATCH --output=/farm_out/%u/%x-%j-%N.out +#SBATCH --error=/farm_out/%u/%x-%j-%N.err +#SBATCH --partition=production +#SBATCH --account=clas12 +#SBATCH --mem-per-cpu=4000 +#SBATCH --gres=disk:1000 +#SBATCH --time=24:00:00 +##SBATCH --mail-user=$USER@jlab.org + +/work/clas12/users/$USER/eic3/largex-eic/macro/dis-5x41/job.sh diff --git a/tutorial/analysis_xqbins.C b/tutorial/analysis_xqbins.C index 3b5f9ac2..309fc36c 100644 --- a/tutorial/analysis_xqbins.C +++ b/tutorial/analysis_xqbins.C @@ -20,7 +20,7 @@ void analysis_xqbins( crossingAngle, outfilePrefix ); - + A->NBINS = 50; // use this to set the number of bins along each axis for each overall bin in e.g. x and Q2 //A->maxEvents = 30000; // use this to limit the number of events A->SetReconMethod("Ele"); // set reconstruction method A->AddFinalState("pipTrack"); // pion final state diff --git a/tutorial/analysis_xqbins_stddevs.C b/tutorial/analysis_xqbins_stddevs.C new file mode 100644 index 00000000..8b986a78 --- /dev/null +++ b/tutorial/analysis_xqbins_stddevs.C @@ -0,0 +1,145 @@ +R__LOAD_LIBRARY(Largex) + +/* run in a grid of (x,Q2) 2D bins + * - various ways to make a grid are demonstrated + * - observe how the resulting histograms differ in each (x,Q2) bin + */ +void analysis_xqbins_stddevs( + TString infiles="datarec/tutorial.config", /* list of input files */ + Double_t eleBeamEn=5, /* electron beam energy [GeV] */ + Double_t ionBeamEn=41, /* ion beam energy [GeV] */ + Double_t crossingAngle=0, /* crossing angle [mrad] */ + TString outfilePrefix="tutorial.xqbins.stddevs." /* output filename prefix*/ +) { + + // setup analysis ======================================== + AnalysisDelphes *A = new AnalysisDelphes( + infiles, + eleBeamEn, + ionBeamEn, + crossingAngle, + outfilePrefix + ); + A->NBINS = 3; // use this to set the number of bins along each axis for each overall bin in e.g. x and Q2 + //A->maxEvents = 30000; // use this to limit the number of events + A->SetReconMethod("Ele"); // set reconstruction method + A->AddFinalState("pipTrack"); // pion final state + //A->AddFinalState("KpTrack"); // kaon final state + //A->AddFinalState("jet"); // jets + + + // define cuts ==================================== + /* - cuts are defined the same way as bins are defined (see below for syntax details and examples); + * the `CutDef` class handles bin definitions + * + * For example, if you want to apply the y<0.95 cut, do: + * A->AddBinScheme("y"); + * A->BinScheme("y")->BuildBin("Max",0.95); + * This makes a single y bin, defined by a maximum of 0.95. + * + * If instead you want to make a cut of 0.01AddBinScheme("y"); + * A->BinScheme("y")->BuildBin("Range",0.01,0.95); + * to define a single y bin. If instead you try to do something like + * A->AddBinScheme("y"); + * A->BinScheme("y")->BuildBin("Min",0.01); + * A->BinScheme("y")->BuildBin("Max",0.95); + * you would actually be defining two y bins, one with the minimum and a + * second with the maximum, which may not be what you want. + * + * You must also be mindful of what other bins you are defining. Suppose you + * want to apply a Q2>10 GeV2 cut, and then do an analysis in bins of Q2. If + * you do + * A->AddBinScheme("q2"); + * A->BinScheme("q2")->BuildBin("Min",10); // your Q2>19 GeV2 cut + * A->BinScheme("q2")->BuildBins( 5, 10, 100, true ); // 5 bins in range 10-100 GeV, equal width log scale + * you are actually defining 6 Q2 bins: the 5 you specify with `BuildBins` + * plus the one you specified with `BuildBin("Min",10)`. In this case, only + * the third line is needed to apply the Q2>10 GeV2 cut, since your binning + * range starts at 10 GeV2. + */ + // here are some common cuts we typically use for SIDIS: + A->AddBinScheme("w"); A->BinScheme("w")->BuildBin("Min",3.0); // W > 3 GeV + A->AddBinScheme("y"); A->BinScheme("y")->BuildBin("Range",0.01,0.95); // 0.01 < y < 0.95 + A->AddBinScheme("z"); A->BinScheme("z")->BuildBin("Range",0.2,0.9); // 0.2 < z < 0.9 + A->AddBinScheme("xF"); A->BinScheme("xF")->BuildBin("Min",0.0); // xF > 0 + A->AddBinScheme("ptLab"); A->BinScheme("ptLab")->BuildBin("Min",0.1); // pT_lab > 0.1 GeV (tracking limit) + + + + // set binning scheme ==================================== + + // first add what bin schemes you want; you don't have to define + // bins for all of them, but you need to add at least the ones + // you plan to use below + A->AddBinScheme("q2"); + A->AddBinScheme("x"); + A->AddBinScheme("pt"); + + // then add the bins that you want + // - tutorial switch statement: change tutorial number to try out + // the different binning implementations + int tutorialNum = 1; + switch(tutorialNum) { + + case 0: + // 1D binning in Q2, equal width in linear scale + // - arguments of BuildBins are: (numBins, min, max, log-scale bool) + // - alternatively: BuildBins(TAxis*, log-scale bool) + // - log-scale is false by default + A->BinScheme("q2")->BuildBins( 3, 1, 100); + break; + + case 1: + // 3x3 grid of (x,Q2) bins, equal width in logarithmic scale + A->BinScheme("q2")->BuildBins( 3, 1, 100, true ); + A->BinScheme("x")->BuildBins( 3, 0.01, 1, true ); + break; + + case 2: + // alternatively: equal width in linear scale + A->BinScheme("q2")->BuildBins( 3, 1, 100 ); + A->BinScheme("x")->BuildBins( 3, 0.01, 1 ); + break; + + case 3: + // custom 2x2 grid (see `CutDef` class for more cut definitions) + // - arguments of BuildBin are (cutType, a, b), where cutType is + // one of the types given in `../src/CutDef.cxx`, and a and b + // depend on which cutType + // - various cutTypes are exemplified here: + A->BinScheme("q2")->BuildBin("Max",10); // Q2 < 10 GeV2 + A->BinScheme("q2")->BuildBin("Min",10); // Q2 > 10 GeV2 + A->BinScheme("x")->BuildBin("Range", 0.05, 0.2 ); // 0.05 < x < 0.2 + A->BinScheme("x")->BuildBin("CenterDelta", 0.5, 0.1 ); // |x-0.5|<0.1 + break; + + case 4: + // overlapping Q2 bins, by specifying various Q2 minima + // - bins are arbitrary and allowed to be overlapping + A->BinScheme("q2")->BuildBin("Min",1); + A->BinScheme("q2")->BuildBin("Min",10); + A->BinScheme("q2")->BuildBin("Min",50); + break; + + case 5: + // 3D binning: 2x2 grid of (x,Q2) bins, for two pT bins + // - you can add more dimensions, but be careful of the curse + // of dimensionality + A->BinScheme("q2")->BuildBins( 2, 1, 100, true ); + A->BinScheme("x")->BuildBins( 2, 0.01, 1, true ); + A->BinScheme("pt")->BuildBin("Range", 0.2, 0.5 ); + A->BinScheme("pt")->BuildBin("Range", 0.5, 0.8 ); + break; + }; + + + + // perform the analysis ================================== + A->Execute(); + + // for reference, here is a print out of HistosDAG + // - it lists each node, together with its inputs and outputs, which + // indicate the connections between the nodes + //A->GetHistosDAG()->PrintBreadth("HistosDAG Nodes"); +}; diff --git a/tutorial/postprocess_xqbins_draw.C b/tutorial/postprocess_xqbins_draw.C index f82e5fa8..20b55a03 100644 --- a/tutorial/postprocess_xqbins_draw.C +++ b/tutorial/postprocess_xqbins_draw.C @@ -49,6 +49,29 @@ void postprocess_xqbins_draw( P->DrawSingle(H,"phiH_RvG","COLZ"); P->DrawSingle(H,"phiS_RvG","COLZ"); */ + + /* + P->DrawSingle(H,"Q2vsX_q2res","COLZ"); + P->DrawSingle(H,"Q2vsX_xres","COLZ"); + P->DrawSingle(H,"Q2vsX_yres","COLZ"); + P->DrawSingle(H,"Q2vsX_zres","COLZ"); + P->DrawSingle(H,"Q2vsX_pTres","COLZ"); + P->DrawSingle(H,"Q2vsX_phiHres","COLZ"); + */ + + P->DrawSingle(H,"Q2vsX_q2resSD","COLZ"); + P->DrawSingle(H,"Q2vsX_xresSD","COLZ"); + P->DrawSingle(H,"Q2vsX_yresSD","COLZ"); + P->DrawSingle(H,"Q2vsX_zresSD","COLZ"); + P->DrawSingle(H,"Q2vsX_pTresSD","COLZ"); + P->DrawSingle(H,"Q2vsX_phiHresSD","COLZ"); + + P->DrawSummary(H,"Q2vsX_q2resSD","COLZ"); + P->DrawSummary(H,"Q2vsX_xresSD","COLZ"); + P->DrawSummary(H,"Q2vsX_yresSD","COLZ"); + P->DrawSummary(H,"Q2vsX_zresSD","COLZ"); + P->DrawSummary(H,"Q2vsX_pTresSD","COLZ"); + P->DrawSummary(H,"Q2vsX_phiHresSD","COLZ"); } ); diff --git a/tutorial/postprocess_xqbins_stddevs_draw.C b/tutorial/postprocess_xqbins_stddevs_draw.C new file mode 100644 index 00000000..dd130570 --- /dev/null +++ b/tutorial/postprocess_xqbins_stddevs_draw.C @@ -0,0 +1,93 @@ +R__LOAD_LIBRARY(Largex) + +// make kinematics coverage plots, such as eta vs. p in bins of (x,Q2) +void postprocess_xqbins_stddevs_draw( + TString infile="out/tutorial.xqbins.stddevs.root" +) { + + // setup postprocessor ======================================== + // - this will read in the Histos objects and your defined + // binning, to construct a HistosDAG object + PostProcessor *P = new PostProcessor(infile); + + // print DAG ================================================== + // - here we print out the HistosDAG object, to show the structure + // - it lists each node, together with its inputs and outputs, which + // indicate the connections between the nodes + // - observe how the structure changes for the different bins + // you specified in the analysis macro + P->Op()->PrintBreadth("HistosDAG Initial Setup"); + + // lambdas ===================================================== + // - we now define what operations we want to happen while + // traversing through the DAG; this example shows how to define + // a payload operator, which will act on every multidimensional + // bin + + // payload: draw a few plots, using PostProcessor::DrawSingle + // - this lambda expression will be executed on every multidimensional + // bin's Histos object + // - capture the pointer to PostProcessor `P` by reference, so we + // have access to it while the lambda executes + // - lambda arguments can include a Histos pointer, and optinally + // a `NodePath` pointer, which gives binning information; here + // we just need the Histos pointer + // - `PostProcessor::Op()` is just an interface (alias) for the + // `HistosDAG` object; see classes `HistosDAG` and its parent `DAG` + // for available methods + P->Op()->Payload( + [&P](Histos *H) { + P->DrawSingle(H,"Q2vsX","COLZ"); + P->DrawSingle(H,"etaVsP","COLZ"); + /* + P->DrawSingle(H,"x_Res",""); + P->DrawSingle(H,"Q2_Res",""); + P->DrawSingle(H,"y_Res",""); + P->DrawSingle(H,"phiH_Res",""); + P->DrawSingle(H,"phiS_Res",""); + P->DrawSingle(H,"x_RvG","COLZ"); + P->DrawSingle(H,"phiH_RvG","COLZ"); + P->DrawSingle(H,"phiS_RvG","COLZ"); + */ + + /* + P->DrawSingle(H,"Q2vsX_q2res","COLZ"); + P->DrawSingle(H,"Q2vsX_xres","COLZ"); + P->DrawSingle(H,"Q2vsX_yres","COLZ"); + P->DrawSingle(H,"Q2vsX_zres","COLZ"); + P->DrawSingle(H,"Q2vsX_pTres","COLZ"); + P->DrawSingle(H,"Q2vsX_phiHres","COLZ"); + */ + + P->DrawSingle(H,"Q2vsX_q2resSD","COLZ"); + P->DrawSingle(H,"Q2vsX_xresSD","COLZ"); + P->DrawSingle(H,"Q2vsX_yresSD","COLZ"); + P->DrawSingle(H,"Q2vsX_zresSD","COLZ"); + P->DrawSingle(H,"Q2vsX_pTresSD","COLZ"); + P->DrawSingle(H,"Q2vsX_phiHresSD","COLZ"); + + P->DrawSummary(H,"Q2vsX_q2resSD","COLZ"); + P->DrawSummary(H,"Q2vsX_xresSD","COLZ"); + P->DrawSummary(H,"Q2vsX_yresSD","COLZ"); + P->DrawSummary(H,"Q2vsX_zresSD","COLZ"); + P->DrawSummary(H,"Q2vsX_pTresSD","COLZ"); + P->DrawSummary(H,"Q2vsX_phiHresSD","COLZ"); + } + ); + + // print DAG ============================ + // - we only added a payload, we did not restructure the DAG, + // therefore this second printout should have the same structure + // as the first + P->Op()->PrintBreadth("HistosDAG Final Setup"); + + // execution =================================================== + // - after you have defined all your operators, run them: + P->Execute(); + + // finish =================================================== + // - PostProcessor::Finish() is NECESSARY to close output streams + // - output file names and directories will be printed to stdout + P->Finish(); + +};