Eli Lichtblau

GITHUB: https://github.com/EliLichtblau/Honr269L/tree/master

       THE GITHUB IS PRIVATE CONTACT ME @elichtbl@umd.edu if you want access with your umd email account

Exact commands to reach current point:

Log onto your umd cluster account:

i.e honr20_somename@siab-1.umd.edu - Use this extension, not hepcms

In a shell:

ssh -Y hon20_someone@siab-1.umd.edu

If you want to mount the cluster so you can use personal editors

Ubuntu - probably works for all debian systems but I'm not sure

sudo apt-get install sshfs

This will mount the cluster to your local directory

That local directory must be empty or it will fail

sshfs something@siab-1.umd.edu:/home/something . 

I'd recommend putting things in folders, actions here are saved to the cluster:

mkdir someDir 

cd someDir

Copy all files from github to directory

Activate root

source ./setup_root

Run H3_DATARQ.C

root -l H3_DATARQ.C

This command will make one of the necessary root files

To make different files

comment out lines 215 and 216 and un-comment 217,218

root -l H3_DATARQ.C

Both root files now exist:

Run makeClean to see plot we care about

root -l makeClean.cc

When cuts start working, run cuts

root -l cuts.cc

Separate Task

All actions assume ubuntu

The cluster is exceedingly slow, especially for file operations, whether this is a function of my wifi, or the school transferring online, I am going to make a local build of root.

As well as download VSCode for auto-completion - will involve changing compile path for c++

First we will download vscode - because its rather simple:

    sudo snap install --classic code   

Install root

https://root.cern.ch/downloading-root - install latest version of root from website

tar -zxvf downloadedFile

cmake - you must navigate to where the required .txt file is

make

make install

Test installation source bin/this.root.sh -assuming shell is bash if shell is tcsh of course use .csh

From here you must change your vscode c++ settings to:

{

"editor.minimap.enabled": true,

"editor.renderControlCharacters": false,

"window.zoomLevel": 0,

"C_Cpp.updateChannel": "Insiders",

"C_Cpp.autocomplete": "Default",

//"C_Cpp.default.compilerPath": "/home/eli/ROOTBUILD/%HOME/usr/local/",

"C_Cpp.default.compilerPath": "/home/eli/root/root-6.20.00",

"C_Cpp.default.includePath": ["/home/eli/ROOTBUILD/%HOME/usr/local/include"],

"C_Cpp.default.intelliSenseMode": "gcc-x64",

"C_Cpp.intelliSenseEngine": "Default",

}

This is just my path

your paths will be different

but you get the idea.

Separate Task

properly cutting the data:

Cutting data from made histograms is a bad idea!

Change the analysis code itself, there are checks already built in:

const double ZMIN=0.0;

const double ZMAX=1500.0;

const double RMAX=735.0;

Did I spend 40 hours trying to reinvent the wheel? Maybe? 

Just change these values to 20., 1370., 695.

done.

           You make a plot which shows the fraction of tritium data with response less than median value of WIMP response vs S1 signal. In this plot, horizontal axis is S1 counts and vertical axis is fraction.

Double_t xbins[nx+1];

Double_t ybins[ny+1];

for(int i = 0;i<nx+1;i++){

xbins[i]=i;

ybins[i]=i;

}

TH2D* mine = (TH2D*)(GGS1Area_phd_logS2AreaOverS1Area_phd_ss->Clone());

for(int x =1;x<GGS1Area_phd_logS2AreaOverS1Area_phd_ss->GetNbinsX(); x++){

for(int y = 1; y<GGS1Area_phd_logS2AreaOverS1Area_phd_ss->GetNbinsY(); y++){

Double_t by = GGS1Area_phd_logS2AreaOverS1Area_phd_ss->GetBinContent(y);

Double_t binContent = GGS1Area_phd_logS2AreaOverS1Area_phd_ss->GetBinContent(x,y);

Double_t guess = GGS1Area_phd_logS2AreaOverS1Area_phd_ss->GetXaxis()->FindBin(x);

if (mine->GetYaxis()->GetBinUpEdge(y)>fFileMean->Eval(x) ){

mine->SetBinContent(x, y, 0);

}

if (binContent != 0){

//cout<<"x "<<x<<" y "<<y<<endl;

cout<<"("<<x<<","<<y<<")"<<(mine->GetBinContent(x,y))<<endl;

}

//mine->Fill(bx,by);

}

}

cout<<mine->GetBinContent(209, 45)<<endl;

mine->GetXaxis()->SetRangeUser(0,200);

mine->GetYaxis()->SetRangeUser(0,5);

//mine->SetMarkerSize(50);

mine->Draw("same");

fFileMean->Draw("same");

       NITTY GRITTY

The work of February 20th through March 9th is best explained by taking the 1500 line function:

and transforming it to do just what we need with this:

 string hname;

    int ix=0;int iy=0; int ipad =0;

    int nx=1;int ny=1;

    TText *t;

    //get wimpl

    f = new TFile("H3_DataRQ_plotsWIMP.root");f->cd();f->ls();

    g = new TFile("H3_DataRQ_plotsNOTWIMP.root");g->cd();g->ls();

    //get means of graphs and style them

    TGraph *fFileMean = getMean(f, "RecoS1Area_phd_S1Bin", "RecoS2Area_phd_S1Bin");

    TGraph *gFileMean = getMean(g, "RecoS1Area_phd_S1Bin", "RecoS2Area_phd_S1Bin");

    fFileMean->SetLineColor(7);

    gFileMean->SetLineColor(6);

    fFileMean->SetLineWidth(4);

    gFileMean->SetLineWidth(4);

    //get graph from f graph

    const char* graph= "S1Area_phd_logS2Area_phd_ss";

    TH2D* S1Area_phd_logS2AreaOverS1Area_phd_ss =(TH2D*) f->Get(graph);  

    TH2D* GGS1Area_phd_logS2AreaOverS1Area_phd_ss = (TH2D*) g->Get(graph); 

    TCanvas* C25 = new TCanvas("C25","C25",1000,1000); C25->cd();

    //set axises

    S1Area_phd_logS2AreaOverS1Area_phd_ss->GetXaxis()->SetRangeUser(0,200);

    S1Area_phd_logS2AreaOverS1Area_phd_ss->GetYaxis()->SetRangeUser(0,5);

    S1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("colz");

    GGS1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("box same");

    gFileMean->Draw("same");

    fFileMean->Draw("same");

    C25->SaveAs("OverlayWithLines.png");

There are function calls but the code itself is very readable, data abstraction even for basic tasks!

The original code contains a 1500 line function...

99% of my time is trying to put things into functions so I can read it,

part of the issue is I'm not great with c++ and pointers get messed up with roots

automatic garbage collection which is a massive headache when putting things in functions

like massive.

Current task is to make cuts that model those done in: Projected WIMP sensitivity of the LUX-ZEPLIN (LZ) dark matter experiment

It is highly unclear for me right now.

All Hope is Lost:

Things I do not know

• Physics

• Histograms

  • Ratio plots

  • error

  • structure of root 

Culmination of Hell

#include "TROOT.h"

#include "TMath.h"

#include "TFile.h"

#include "TTree.h"

#include "TObject.h"

#include "TH1D.h"

#include "TH2D.h"

#include "TProfile.h"

#include "TF1.h"

#include <TGraph.h>

#include <TGraphErrors.h>

#include <TStyle.h>

#include <TCanvas.h>

#include <THStack.h>

#include <TMinuit.h>

#include <TLatex.h>

#include <iostream>

#include <fstream>

#include <string>

#include <iomanip>

#include <cmath>

#include <vector>

void SetHistStyle(TH1D* h){

}

TFile *wimp;

TFile *notWimp;

int NS1BINS=20;

std::vector<TPad*> make_pads(TCanvas* c, int nx, int ny){

c->cd();

std::vector<TPad*> PadList;

double dx=1.0/float(nx);

double dy=1.0/float(ny);

for(int ix=0;ix<nx;ix++){

for(int iy=0;iy<ny;iy++){

double x1=ix*dx;

double x2=(ix+1)*dx;

double y1=iy*dy;

double y2=(iy+1)*dy;

TPad* pad = new TPad("c","c",x1,y1,x2,y2);

PadList.push_back(pad);

}

}

return PadList;

}

std::vector<TH1D*> GetTH1DVector(int NBINS,string hname1, TFile *f){

std::vector<TH1D*> TH1DVector;

string hname;

for(int i=0;i<NBINS;i++){

hname = hname1+to_string(i);

TH1DVector.push_back((TH1D*) f->Get(hname.c_str()));

}

return TH1DVector;

}

std::vector<TH2D*> GetTH2DVector(int NBINS, string hname1, TFile *f){

std::vector<TH2D*> TH2DVector;

string hname;

for(int i = 0; i < NBINS; i++){

hname = hname1 + to_string(i);

TH2DVector.push_back((TH2D*) f->Get(hname.c_str()));

}

return TH2DVector;

}

std::vector<TProfile*> GetTProfileVector(int NBINS,string hname1, TFile *f){

std::vector<TProfile*> TProfileVector;

string hname;

for(int i=0;i<NBINS;i++){

hname = hname1+to_string(i);

TProfileVector.push_back((TProfile*) f->Get(hname.c_str()));

}

return TProfileVector;

}

void effX() {

auto effErf = [](double* x, double* p) {

return (TMath::Erf((x[0] - p[0]) / p[1]) + 1) / 2. * p[2];

};

TF1* myErf = new TF1("myErf", effErf, 0., 10.);

myErf->SetParameter(0, 5.);

myErf->SetParameter(1, 5.);

myErf->SetParameter(2, 1.);

// eff->Fit(myErf);

}

double FermiFunction(double T){

double masse = 511; // mass of electron

double Zp = 2.0; // charge of daughter element i.e. He

double alpha=1.0/137.; // The alpha_em

double beta=sqrt(2*T/masse); // velocity

double xtmp = 2.0*M_PI*alpha*Zp/beta;

double denom = (1.0-exp(-xtmp));

double fermi = xtmp/denom;

return fermi;

}

double TritiumDecaySpectrum(double T){

double Q=18.590;

if(T>Q) return 0;

if(T<=0) return 0;

double m_nue=0.;

double m_electron=511.0;

// double Fermi=1.0;

double Norm=1.0e-5/1.00121;

double C=1.0e-3;

double a1 = sqrt(pow(T,2.)+2*T*m_electron);

double a2 = T+m_electron;

double a3 = Q-T;

double a4 = sqrt(pow(a3,2.)-pow(m_nue,2.));

double weight = Norm*C*a1*a2*a3*a4*FermiFunction(T);

return weight;

}

Double_t H3Spectrum(Double_t *x,Double_t *par){

double T= x[0];

Double_t fitval = par[0]*TritiumDecaySpectrum(T);

return fitval;

}

TGraph *getMean(TFile *f, string xVectorHist, string yVectorHist) {

TGraph *mean;

//get x and y vectors

std::vector<TH1D*> RecoS1AreaVec = GetTH1DVector(NS1BINS,xVectorHist, f);

std::vector<TH1D*> RecoS2Area_phdOverS1AreaVec = GetTH1DVector(16,yVectorHist, f);

//get x entries

int npoints = RecoS2Area_phdOverS1AreaVec.size();

//this is literally the same as the class reference, still dont understand

const int nq=100;

double xq[nq];

for(int i=0;i<nq;i++){

xq[i]= double(i+1)/nq;

}

//S1 mean values over npoints

double S1Mean[npoints];

double yq[npoints][nq];

for(int i=0;i<npoints;i++){

//set mean value to reconstructed values

S1Mean[i]= RecoS1AreaVec[i]->GetMean();

//sets yq to quantile values

RecoS2Area_phdOverS1AreaVec[i]->GetQuantiles(nq,yq[i],xq);

}

double logQ50[npoints];

for(int i=0;i<npoints;i++){

logQ50[i] = log10(yq[i][50]);

}

mean = new TGraph(npoints,S1Mean,logQ50);

return mean;

}

double* getMedian(TFile *f, string xVectorHist, string yVectorHist, double* S1Mean) {

TGraph *mean;

//get x and y vectors

std::vector<TH1D*> RecoS1AreaVec = GetTH1DVector(NS1BINS,xVectorHist, f);

std::vector<TH1D*> RecoS2Area_phdOverS1AreaVec = GetTH1DVector(16,yVectorHist, f);

//get x entries

int npoints = RecoS2Area_phdOverS1AreaVec.size();

//this is literally the same as the class reference, still dont understand

const int nq=100;

double xq[nq];

for(int i=0;i<nq;i++){

xq[i]= double(i+1)/nq;

}

//S1 mean values over npoints

double yq[npoints][nq];

for(int i=0;i<npoints;i++){

//set mean value to reconstructed values

S1Mean[i]= RecoS1AreaVec[i]->GetMean();

//sets yq to quantile values

RecoS2Area_phdOverS1AreaVec[i]->GetQuantiles(nq,yq[i],xq);

}

return S1Mean;

}

void s1Xevents_s2XEvents(TFile *f, TCanvas* c, int colorS1, int colorS2){

c->cd();

auto s1_area_phd = (TH1D*) f->Get("S1Area_phd_ss");

auto s2_area_phd = (TH1D*) f->Get("S2Area_phd_ss");

cout<<s1_area_phd<<endl;

cout<<s2_area_phd<<endl;

cout<<"\n\n\n";

s1_area_phd->GetXaxis()->SetRangeUser(0,5e5);

s1_area_phd->GetYaxis()->SetRangeUser(0,2000);

s2_area_phd->GetXaxis()->SetRangeUser(0,5e5);

s2_area_phd->GetYaxis()->SetRangeUser(0,2000);

s1_area_phd->SetLineColor(colorS1);

s2_area_phd->SetLineColor(colorS2);

s2_area_phd->SetLineWidth(2);

s1_area_phd->SetLineWidth(2);

//cout<<"t3 "<<t3<<endl;

//t3->DrawClone("same");

//temp->DrawClone("same");

s1_area_phd->Draw("same");

s2_area_phd->Draw("same");

}

void ratio_plot(){

string hname;

int ix=0;int iy=0; int ipad =0;

int nx=1;int ny=1;

TText *t;

//get wimp and not wimp data

wimp = new TFile("H3_DataRQ_plotsWIMPCut.root");wimp->cd();//wimp->ls();

notWimp = new TFile("H3_DataRQ_plotsNOTWIMPCut.root");notWimp->cd();//notWimp->ls();

TH1D* s1_area_phd = (TH1D*) wimp->Get("S1Area_phd_ss");

TH1D* s1_area_phd_n = (TH1D*) notWimp->Get("S1Area_phd_ss");

TH1D* s2_area_phd = (TH1D*) wimp->Get("S2Area_phd_ss");

//get graph from f graph

const char* graph= "S1Area_phd_logS2Area_phd_ss";

TH2D* S1Area_phd_logS2AreaOverS1Area_phd_ss =(TH2D*) wimp->Get(graph);

TH2D* GGS1Area_phd_logS2AreaOverS1Area_phd_ss = (TH2D*) notWimp->Get(graph);

//set axises

S1Area_phd_logS2AreaOverS1Area_phd_ss->GetXaxis()->SetRangeUser(0,200);

S1Area_phd_logS2AreaOverS1Area_phd_ss->GetYaxis()->SetRangeUser(0,5);

TCanvas* C3 = new TCanvas("C3","C3",1000,1000); C3->cd();

/*medianWimp response across s1bins*/

std::vector<TH1D*> RecoS1AreaVec = GetTH1DVector(NS1BINS,"RecoS1Area_phd_S1Bin", wimp);

std::vector<TH1D*> RecoS2Area_phdOverS1AreaVec = GetTH1DVector(16,"RecoS2Area_phd_S1Bin", wimp);

//get x entries

int npoints = RecoS2Area_phdOverS1AreaVec.size();

//this is literally the same as the class reference, still dont understand

const int nq=100;

double xq[nq];

for(int i=0;i<nq;i++){

xq[i]= double(i+1)/nq;

}

//S1 mean values over npoints

double medianWimp[npoints];

double amountBelow[npoints];

double amountTotal[npoints];

for(int i=0;i<npoints;i++){

medianWimp[i]= RecoS1AreaVec[i]->GetMean();

for(int x = 0; x < RecoS1AreaVec[i]->GetXaxis()->GetNbins();x++){

int binval = RecoS1AreaVec[i]->GetBinCenter(x);

int binfreq = RecoS1AreaVec[i]->GetBinContent(x);

if (binval < medianWimp[i]){

amountBelow[i] += binfreq;

}

amountTotal[i] += binfreq;

}

cout<<amountBelow[i]/amountTotal[i]<<endl;

}

//for s1 bin calculate number of H3 response < median[s1bin]

Double_t ratioplotArr[npoints+1];

for(int i =0;i<npoints+1;i++){ratioplotArr[i]=i;}

TH1D* ratioPlot = new TH1D("ratioplot", "ratioplot ", npoints, ratioplotArr);

for(int i =0;i<npoints+1;i++){

ratioPlot->SetBinContent(i, amountBelow[i]/amountTotal[i]);

ratioPlot->SetBinError(i, (amountBelow[i]/amountTotal[i])/10);

}

//R[s1bin] = H3pass[s1bin]/H3Total[s1bin

//fill histogram setbincontent(R[s1bin], s1bin)... setbinerror

ratioPlot->Draw("same");

//S1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("box same");

//GGS1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("same");

C3->SaveAs("working onit.png");

}

string hname;

int ix=0;int iy=0; int ipad =0;

int nx=1;int ny=1;

TText *t;

//get wimp and not wimp data

wimp = new TFile("H3_DataRQ_plotsWIMPCut.root");wimp->cd();//wimp->ls();

notWimp = new TFile("H3_DataRQ_plotsNOTWIMPCut.root");notWimp->cd();//notWimp->ls();

//GetTH1DVector(int NBINS,string hname1, TFile *f)

//std::vector<TH1D*> wimpRawS1 = GetTH1DVector(500, "RawS1Photons_ss", wimp);

//std::vector<TH1D*> notWimpRawS1 = GetTH1DVector(NS1BINS, "RawS2Photons_ss", notWimp);

TCanvas* C1 = new TCanvas("C1","C1",1000,1000); C1->cd();

nx=1;ny=2; std::vector<TPad*> C1Pads = make_pads(C1,nx,ny);

TH1D* s1_area_phd = (TH1D*) wimp->Get("S1Area_phd_ss");

TH1D* s1_area_phd_n = (TH1D*) notWimp->Get("S1Area_phd_ss");

TH1D* s2_area_phd = (TH1D*) wimp->Get("S2Area_phd_ss");

//s1Xevents_s2XEvents(wimp,C1, 2,3);

//c->cd();

cout<<s1_area_phd<<endl;

cout<<s2_area_phd<<endl;

s1_area_phd->SetLineColor(2);

s2_area_phd->SetLineColor(3);

s2_area_phd->SetLineWidth(2);

s1_area_phd->SetLineWidth(2);

//s2_area_phd->GetXaxis()->SetRangeUser(0,5e5);

//s2_area_phd->GetYaxis()->SetRangeUser(0,2000);

//cout<<"t3 "<<t3<<endl;

//t3->DrawClone("same");

//temp->DrawClone("same");

//s1_area_phd->Draw();

//idea->Draw();

ix=0;iy=0; ipad = nx*ix+iy; C1->cd();C1Pads[ipad]->Draw(); C1Pads[ipad]->cd(); s1_area_phd->Draw("same");

ix=0;iy=0; ipad = nx*ix+iy; C1->cd();C1Pads[ipad]->Draw(); C1Pads[ipad]->cd(); s1_area_phd_n->Draw("same");

ix=0;iy=1; ipad = nx*ix+iy; C1->cd();C1Pads[ipad]->Draw(); C1Pads[ipad]->cd(); s2_area_phd->Draw("same");

//s1_area_phd->Draw("same");

//s2_area_phd->Draw("same");

//cout<<s1->getColor()<<endl;

//s1Xevents_s2XEvents(notWimp,C1,4,5);

//C1->SaveAs("RawOverlay.png");

//get means of graphs and style them

TGraph *fFileMean = getMean(wimp, "RecoS1Area_phd_S1Bin", "RecoS2Area_phd_S1Bin");

TGraph *gFileMean = getMean(notWimp, "RecoS1Area_phd_S1Bin", "RecoS2Area_phd_S1Bin");

fFileMean->SetLineColor(7);

gFileMean->SetLineColor(6);

fFileMean->SetLineWidth(1);

gFileMean->SetLineWidth(4);

//get graph from f graph

const char* graph= "S1Area_phd_logS2Area_phd_ss";

TH2D* S1Area_phd_logS2AreaOverS1Area_phd_ss =(TH2D*) wimp->Get(graph);

TH2D* GGS1Area_phd_logS2AreaOverS1Area_phd_ss = (TH2D*) notWimp->Get(graph);

TCanvas* C25 = new TCanvas("C25","C25",1000,1000); C25->cd();

//set axises

S1Area_phd_logS2AreaOverS1Area_phd_ss->GetXaxis()->SetRangeUser(0,200);

S1Area_phd_logS2AreaOverS1Area_phd_ss->GetYaxis()->SetRangeUser(0,5);

S1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("colz");

GGS1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("box same");

gFileMean->Draw("same");

fFileMean->Draw("same");

C25->SaveAs("OverlayWithLines.png");

TCanvas* C3 = new TCanvas("C3","C3",1000,1000); C3->cd();

TH2D* mine = (TH2D*)(S1Area_phd_logS2AreaOverS1Area_phd_ss->Clone());

/*for(int x =0;x<mine->GetNbinsX(); x++){

for(int y = 0; y<mine->GetNbinsY(); y++){

if (mine->GetYaxis()->GetBinCenter(y)>fFileMean->Eval(x) ){

mine->SetBinContent(x, y, 0);

}

}

}*/

int sizeX = mine->GetNbinsX();

double conditonal[sizeX];

double total[sizeX];

double ratio[sizeX];

for(int x =0;x<sizeX; x++){

Int_t xBinTotal = 0;

Int_t conditionXBin = 0;

for(int y = 0; y<mine->GetNbinsY(); y++){

xBinTotal += mine->GetBinContent(x,y);

if (mine->GetYaxis()->GetBinCenter(y)>fFileMean->Eval(x) ){

//mine->SetBinContent(x, y, 0);

conditionXBin += mine->GetBinContent(x,y);

}

}

total[x] = xBinTotal;

if (xBinTotal == 0){

ratio[x] = 0.;

}else{

ratio[x] = ((Float_t)conditionXBin)/xBinTotal;

}

}

Double_t ratHist[sizeX+1];

for(int i =0;i<=sizeX;i++){ratHist[i]=i;}

TH1D* hist_ratio = new TH1D("histRatio", "histRatio", sizeX, ratHist);

TH1D* hist_S1count = new TH1D("histRatio", "histRatio", sizeX, ratHist);

for(int i =0;i<=sizeX;i++){

hist_ratio->SetBinContent(i,ratio[i]);

hist_S1count->SetBinContent(i,total[i]);

}

hist_ratio->GetYaxis()->SetTitle("Fraction");

hist_S1count->GetYaxis()->SetTitle("s1 count");

TRatioPlot* rp = new TRatioPlot(hist_ratio, hist_S1count);

rp->Draw();

//total_over_conditional->Draw();

//cout<<mine->GetBinContent(209, 45)<<endl;

mine->GetXaxis()->SetRangeUser(0,200);

mine->GetYaxis()->SetRangeUser(0,5);

mine->SetMarkerSize(50);

mine->Draw("same");

fFileMean->Draw("same");

//S1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("same colz");

C3->SaveAs("working onit.png");

string hname;

int ix=0;int iy=0; int ipad =0;

int nx=1;int ny=1;

TText *t;

//get wimp and not wimp data

DataFile = new TFile("BG_DataRQ_plots.root");DataFile->cd();//wimp->ls();

//GetTH1DVector(int NBINS,string hname1, TFile *f)

//std::vector<TH1D*> wimpRawS1 = GetTH1DVector(500, "RawS1Photons_ss", wimp);

//std::vector<TH1D*> notWimpRawS1 = GetTH1DVector(NS1BINS, "RawS2Photons_ss", notWimp);

TCanvas* C1 = new TCanvas("C1","C1",1000,1000); C1->cd();

nx=1;ny=2; std::vector<TPad*> C1Pads = make_pads(C1,nx,ny);

TH1D* correctedX = (TH1D*) DataFile->Get("correctedX_cm_0");

TH1D* correctedY = (TH1D*)DataFile->Get("correctedY_cm_0");

TH1D* driftTime = (TH1D*)DataFile->Get("drift_time");

//style bullshit

correctedX->SetLineColor(8);

correctedY->SetLineColor(11);

//correctedX->Draw("same");

//correctedY->Draw("same");

int nbins = correctedY->GetNbinsX();

cout<<correctedX->GetNbinsX()<<endl;

cout<<correctedY->GetNbinsX()<<endl;

//s2_area_phd->GetXaxis()->SetRangeUser(0,5e5);

//s2_area_phd->GetYaxis()->SetRangeUser(0,2000);

TH1D theboy = (*correctedY)*(*correctedY) + (*correctedX)*(*correctedX);

TH1D* aspointer = (TH1D*)theboy.Clone();

//aspointer->Draw("same");

std::map<double, double> cmXMap;

std::map<double, double> cmYMap;

std::map<double, double> cmXsquare;

std::map<double, double> cmYsquare;

for(int i =0; i<nbins;i++){

double xval = correctedX->GetXaxis()->GetBinCenter(i);

double xfreq = correctedX->GetBinContent(i);

cmXMap[xval] = xfreq;

double yval = correctedY->GetXaxis()->GetBinCenter(i);

double yfreq = correctedY->GetBinContent(i);

cmYMap[yval] = yfreq;

cmXsquare[xval*xval] = xfreq;

cmYsquare[yval*yval] = yfreq;

}

std::map<double, double> cmRsquared;

for(auto& yval : cmYsquare){

//this should get frequency for some yval

//double tmp = cmRsquared[yval.first];

cmRsquared[yval.first] = cmYsquare[yval.first]+cmXsquare[yval.first];

/*if (tmp == NULL){

cmRsquared[yval.first] = cmYsquare[yval.first]+cmXsquare[yval.first];

} else{

cmRsquared[yval.first] = tmp+ cmYsquare[yval.first]+cmXsquare[yval.first];

}*/

}

cout<<cmRsquared.size()<<endl;

double rSquared[cmRsquared.size()];

int index = 0;

int cmRmin;

int cmRmax;

for(auto& val : cmRsquared){

if (index==0){

cmRmin = val.first;

}

rSquared[index] = (val.first);

cmRmax = val.first;

index++;

}

TH1D* myR = new TH1D("rsquare", "tile", cmRsquared.size()-1, rSquared);

myR->SetLineColor(1);

cout<<cmRmin<< " "<<cmRmax<<endl;

index = 0;

for(auto& val : cmRsquared){

//myR->Fill(Double_t(val.first), Double_t(val.second));

myR->Fill(Double_t((val.first)),Double_t(val.second));

index++;

}

myR->GetXaxis()->SetRangeUser(0,10000);

myR->SetLineStyle(1);

myR->SetLineColor(2);

myR->UseCurrentStyle();

myR->Draw("same");

C1->SaveAs("AAAAAAAHFIOEAJFNIOEANFEA.png");

string hname;

int ix=0;int iy=0; int ipad =0;

int nx=1;int ny=1;

TText *t;

//get wimpl

f = new TFile("H3_DataRQ_plotsWIMP.root");f->cd();f->ls();

g = new TFile("H3_DataRQ_plotsNOTWIMP.root");g->cd();g->ls();

//get means of graphs and style them

TGraph *fFileMean = getMean(f, "RecoS1Area_phd_S1Bin", "RecoS2Area_phd_S1Bin");

TGraph *gFileMean = getMean(g, "RecoS1Area_phd_S1Bin", "RecoS2Area_phd_S1Bin");

fFileMean->SetLineColor(7);

gFileMean->SetLineColor(6);

fFileMean->SetLineWidth(4);

gFileMean->SetLineWidth(4);

//get graph from f graph

const char* graph= "S1Area_phd_logS2Area_phd_ss";

TH2D* S1Area_phd_logS2AreaOverS1Area_phd_ss =(TH2D*) f->Get(graph);

TH2D* GGS1Area_phd_logS2AreaOverS1Area_phd_ss = (TH2D*) g->Get(graph);

TCanvas* C25 = new TCanvas("C25","C25",1000,1000); C25->cd();

//set axises

S1Area_phd_logS2AreaOverS1Area_phd_ss->GetXaxis()->SetRangeUser(0,200);

S1Area_phd_logS2AreaOverS1Area_phd_ss->GetYaxis()->SetRangeUser(0,5);

S1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("colz");

GGS1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("box same");

gFileMean->Draw("same");

fFileMean->Draw("same");

C25->SaveAs("OverlayWithLines.png");

//clean up

delete C25;

delete S1Area_phd_logS2AreaOverS1Area_phd_ss;

delete GGS1Area_phd_logS2AreaOverS1Area_phd_ss;

delete gFileMean;

delete fFileMean;

delete f;

delete g;

string hname;

int ix=0;int iy=0; int ipad =0;

int nx=1;int ny=1;

TText *t;

gStyle->SetOptStat(kTRUE);

gStyle->SetOptStat(11111111);

gStyle->SetOptTitle(kTRUE);

//f = new TFile("H3_DataRQ_plots.root");f->cd();f->ls();

f = new TFile("H3_DataRQ_plots.root");f->cd();f->ls();

//wimp = new TFile("H3_DataRQ_plotsWIMP.root");wimp->cd();wimp->ls();

//notWimp = new TFile("H3_DataRQ_plotsNOTWIMP.root");notWimp->cd();notWimp->ls();

TH2D* ParentPositionRSqrZAll = (TH2D*) f->Get("ParentPositionRSqrZAll");

TH2D* ParentPositionRSqrZ_All = (TH2D*) f->Get("ParentPositionRSqrZ_All");

TH2D* ParentPositionRSqrZ_ss = (TH2D*) f->Get("ParentPositionRSqrZ_ss");

TH2D* ParentPositionRSqrZ_ms = (TH2D*) f->Get("ParentPositionRSqrZ_ms");

TH2D* ParentPositionRSqrZ_Kr83m = (TH2D*) f->Get("ParentPositionRSqrZ_Kr83m");

TH2D* ParentPositionRSqrZ_Pileup = (TH2D*) f->Get("ParentPositionRSqrZ_Pileup");

TH2D* ParentPositionRSqrZ_Other = (TH2D*) f->Get("ParentPositionRSqrZ_Other");

/*

fjeaoifjoeajnfoieiahfoeafbneaoihfcioea

*/

hname = "RecoS1Area_phd_S1Bin";

std::vector<TH1D*> RecoS1AreaVec = GetTH1DVector(NS1BINS,hname);

std::vector<TH1D*> RecEnergy_keVVec = GetTH1DVector(NRecoEnergyBins,hname);

int npoints = RecEnergy_keVVec.size()-1;

hname = "RecoS2Area_phdOverS1Area_phd_S1Bin";

std::vector<TH1D*> RecoS2Area_phdOverS1AreaVec = GetTH1DVector(16,hname);

npoints = RecoS2Area_phdOverS1AreaVec.size();

double S1Mean[npoints];

double logS2OverS1Mean[npoints];

double logS2OverS1StdDev[npoints];

std::vector<double> S1Mean0;

const int nq=100;

double xq[nq];

double yq[npoints][nq];

for(int i=0;i<npoints;i++){

S1Mean[i]= RecoS1AreaVec[i]->GetMean();

//S1Mean0.push_back(RecoS1AreaVec[i]->GetMean());

logS2OverS1Mean[i] = log10(RecoS2Area_phdOverS1AreaVec[i]->GetMean());

logS2OverS1StdDev[i] = log10(RecoS2Area_phdOverS1AreaVec[i]->GetStdDev());

RecoS2Area_phdOverS1AreaVec[i]->GetQuantiles(nq,yq[i],xq);

}

double ex[20];

double ey[20];

double Q02[npoints];

double Q05[npoints];

double Q10[npoints];

double Q32[npoints];

double Q50[npoints];

double Q68[npoints];

double Q90[npoints];

double Q95[npoints];

double Q98[npoints];

double logQ02[npoints];

double logQ05[npoints];

double logQ10[npoints];

double logQ32[npoints];

double logQ50[npoints];

double logQ68[npoints];

double logQ90[npoints];

double logQ95[npoints];

double logQ98[npoints];

for(int i=0;i<npoints;i++){

int q;

q=2; Q02[i] = yq[i][q];

q=5; Q05[i] = yq[i][q];

q=10; Q10[i] = yq[i][q];

q=32; Q32[i] = yq[i][q];

q=50; Q50[i] = yq[i][q];

q=68; Q68[i] = yq[i][q];

q=90; Q90[i] = yq[i][q];

q=95; Q95[i] = yq[i][q];

q=98; Q98[i] = yq[i][q];

}

for(int i=0;i<npoints;i++){

int q;

q=2; logQ02[i] = log10(yq[i][q]);

q=5; logQ05[i] = log10(yq[i][q]);

q=10; logQ10[i] = log10(yq[i][q]);

q=32; logQ32[i] = log10(yq[i][q]);

q=50; logQ50[i] = log10(yq[i][q]);

q=68; logQ68[i] = log10(yq[i][q]);

q=90; logQ90[i] = log10(yq[i][q]);

q=95; logQ95[i] = log10(yq[i][q]);

q=98; logQ98[i] = log10(yq[i][q]);

}

TGraphErrors* logS2OverS1Graph = new TGraphErrors(npoints,S1Mean,logS2OverS1Mean,ex,ey);

TH2D* S1Area_phd_S2Area_phd_ss = (TH2D*) f->Get("S1Area_phd_S2Area_phd_ss");

TH2D* S1Area_phd_logS2AreaOverS1Area_phd_ss = (TH2D*) f->Get("S1Area_phd_logS2AreaOverS1Area_phd_ss");

TGraph *gr02;

TGraph *gr05;

TGraph *gr10;

TGraph *gr32;

TGraph *gr50;

TGraph *gr68;

TGraph *gr90;

TGraph *gr95;

TGraph *gr98;

gr02 = new TGraph(npoints,S1Mean,logQ02);

gr05 = new TGraph(npoints,S1Mean,logQ05);

gr10 = new TGraph(npoints,S1Mean,logQ10);

gr32 = new TGraph(npoints,S1Mean,logQ32);

gr50 = new TGraph(npoints,S1Mean,logQ50);

gr68 = new TGraph(npoints,S1Mean,logQ68);

gr90 = new TGraph(npoints,S1Mean,logQ90);

gr95 = new TGraph(npoints,S1Mean,logQ95);

gr98 = new TGraph(npoints,S1Mean,logQ98);

gr02->SetLineWidth(4);

gr05->SetLineWidth(4);

gr10->SetLineWidth(4);

gr32->SetLineWidth(4);

gr50->SetLineWidth(4);

gr68->SetLineWidth(4);

gr90->SetLineWidth(4);

gr95->SetLineWidth(4);

gr98->SetLineWidth(4);

gr02->SetLineColor(7);

gr05->SetLineColor(6);

gr10->SetLineColor(4);

gr32->SetLineColor(3);

gr50->SetLineColor(2);

gr68->SetLineColor(3);

gr90->SetLineColor(4);

gr95->SetLineColor(6);

gr98->SetLineColor(7);

gr98->SetMaximum(4.5);

gr98->SetMinimum(1.5);

logS2OverS1Graph->SetMaximum(4.0);

logS2OverS1Graph->SetMinimum(2.0);

/*

fjeaoifjoeajnfoieiahfoeafbneaoihfcioea

*/

TCanvas* C25 = new TCanvas("C25","C25",1000,1000); C25->cd();

nx=1;ny=1; std::vector<TPad*> C25Pads = make_pads(C25,nx,ny);

S1Area_phd_S2Area_phd_ss->GetXaxis()->SetRangeUser(0,200);

S1Area_phd_S2Area_phd_ss->GetYaxis()->SetRangeUser(0,2e5);

S1Area_phd_logS2AreaOverS1Area_phd_ss->GetXaxis()->SetRangeUser(0,200);

S1Area_phd_logS2AreaOverS1Area_phd_ss->GetYaxis()->SetRangeUser(0,10);

S1Area_phd_logS2AreaOverS1Area_phd_ss->GetXaxis()->SetRangeUser(0,200);

S1Area_phd_logS2AreaOverS1Area_phd_ss->GetYaxis()->SetRangeUser(2.0,4.0);

//S1Area_phd_logS2AreaOverS1Area_phd_ss->GetYaxis()->SetRangeUser(1.0, 5.0);

// ix=0;iy=1; ipad = ny*ix+iy; C25->cd();C25Pads[ipad]->Draw(); C25Pads[ipad]->cd(); S1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("colz");

ix=0;iy=0; ipad = ny*ix+iy; C25->cd();C25Pads[ipad]->Draw(); C25Pads[ipad]->cd(); S1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("colz");

t = new TText(5.0,0.2,"sample text");

t->SetTextAlign(11);

t->SetTextColor(kRed+2);

t->SetTextFont(43);

t->SetTextSize(40);

t->SetTextAngle(0);

// t->Draw();

logS2OverS1Graph->SetMarkerSize(1.0);

logS2OverS1Graph->Draw("P same");

gr98->Draw("same");

gr95->Draw("same");

gr90->Draw("same");

gr68->Draw("same");

gr50->Draw("same");

gr32->Draw("same");

gr10->Draw("same");

gr05->Draw("same");

gr02->Draw("same");

gr98->Draw("same");

C25->SaveAs("NOTWIMP.png");

Very little was accomplished, much toil was indulged, someone once said when c++ is your hammer every problem looks like a thumb - I like c++, but when root is your hammer, every nail is a pressure sensitive landmine that you have no choice but to hit no matter how much you don't want to.

In the past 3 days I have dropped 6lbs fixing code, except I haven't been wrong, my group sent me the wrong data....

The code Anwar Sent:

#include "TROOT.h"

#include "TMath.h"

#include "TFile.h"

#include "TTree.h"

#include "TObject.h"

#include "TH1D.h"

#include "TH2D.h"

#include "TProfile.h"

#include "TF1.h"

#include <TGraph.h>

#include <TGraphErrors.h>

#include <TStyle.h>

#include <TCanvas.h>

#include <THStack.h>

#include <TMinuit.h>

#include <TLatex.h>

#include <iostream>

#include <fstream>

#include <string>

#include <iomanip>

#include <cmath>

#include <vector>

void SetHistStyle(TH1D* h){

}

TFile *f;

TGraphErrors* H3Acceptance;

TGraphErrors* GrThresholdError;

TGraph* GrThreshold;

std::vector<double> diffRecoParentEnergyParent;

std::vector<double> diffRecoParentEnergyNorm;

std::vector<double> diffRecoParentEnergyMean;

std::vector<double> diffRecoParentEnergyStdDev;

int NEBINS=20;

int NS1BINS=20;

int NRecoEnergyBins=20;

TH1D* tspectrum;

TH1D* tspectrum_all;

TH1D* tspectrum_ss;

//TH1D* energyParentAll;

TH1D* energyParent_All;

TH1D* energyParent_ss;

TH1D* energyRecBeforeFit_keV;

//TH1D* energyParentBeforeFit_keV;

std::vector<TPad*> make_pads(TCanvas* c, int nx, int ny){

c->cd();

std::vector<TPad*> PadList;

double dx=1.0/float(nx);

double dy=1.0/float(ny);

for(int ix=0;ix<nx;ix++){

for(int iy=0;iy<ny;iy++){

double x1=ix*dx;

double x2=(ix+1)*dx;

double y1=iy*dy;

double y2=(iy+1)*dy;

TPad* pad = new TPad("c","c",x1,y1,x2,y2);

PadList.push_back(pad);

}

}

return PadList;

}

std::vector<TH1D*> GetTH1DVector(int NBINS,string hname1){

std::vector<TH1D*> TH1DVector;

string hname;

for(int i=0;i<NBINS;i++){

hname = hname1+to_string(i);

TH1DVector.push_back((TH1D*) f->Get(hname.c_str()));

}

return TH1DVector;

}

std::vector<TProfile*> GetTProfileVector(int NBINS,string hname1){

std::vector<TProfile*> TProfileVector;

string hname;

for(int i=0;i<NBINS;i++){

hname = hname1+to_string(i);

TProfileVector.push_back((TProfile*) f->Get(hname.c_str()));

}

return TProfileVector;

}

void effX() {

auto effErf = [](double* x, double* p) {

return (TMath::Erf((x[0] - p[0]) / p[1]) + 1) / 2. * p[2];

};

TF1* myErf = new TF1("myErf", effErf, 0., 10.);

myErf->SetParameter(0, 5.);

myErf->SetParameter(1, 5.);

myErf->SetParameter(2, 1.);

// eff->Fit(myErf);

}

double FermiFunction(double T){

double masse = 511; // mass of electron

double Zp = 2.0; // charge of daughter element i.e. He

double alpha=1.0/137.; // The alpha_em

double beta=sqrt(2*T/masse); // velocity

double xtmp = 2.0*M_PI*alpha*Zp/beta;

double denom = (1.0-exp(-xtmp));

double fermi = xtmp/denom;

return fermi;

}

/*

Minimizer is Linear

Chi2 = 249.85

NDf = 179

p0 = 3.19291 +/- 0.00741538

p1 = 0.0588292 +/- 0.00138176

p2 = -0.00219583 +/- 8.59119e-05

p3 = 4.72473e-05 +/- 2.43349e-06

p4 = -5.74813e-07 +/- 3.5487e-08

p5 = 3.90672e-09 +/- 2.75161e-10

p6 = -1.37867e-11 +/- 1.07518e-12

p7 = 1.96281e-14 +/- 1.66177e-15

*/

double GetMedianWimp_old(double x) {

double p0 = 3.19291;

double p1 = 0.0588292;

double p2 = -0.00219583;

double p3 = 4.72473e-05;

double p4 = -5.74813e-07;

double p5 = 3.90672e-09;

double p6 = -1.37867e-11;

double p7 = 1.96281e-14;

double y=p0+p1*pow(x,1)+p2*pow(x,2)+p3*pow(x,3)+p4*pow(x,4)+p5*pow(x,5)+p6*pow(x,6)+p7*pow(x,7);

return y;

}

/*

inimizer is Linear

Chi2 = 1471.05

NDf = 351

p0 = 3.12198 +/- 0.0109153

p1 = 0.0818549 +/- 0.00290273

p2 = -0.00454078 +/- 0.000269032

p3 = 0.00015603 +/- 1.195e-05

p4 = -3.26929e-06 +/- 2.9261e-07

p5 = 4.23779e-08 +/- 4.2257e-09

p6 = -3.39294e-10 +/- 3.6816e-11

p7 = 1.62543e-12 +/- 1.89769e-13

p8 = -4.25605e-15 +/- 5.31928e-16

p9 = 4.67382e-18 +/- 6.23991e-19

*/

double GetMedianWimp(double x) {

double p0 = 3.12198;

double p1 = 0.0818549;

double p2 = -0.00454078;

double p3 = 0.00015603;

double p4 = -3.26929e-06;

double p5 = 4.23779e-08;

double p6 = -3.39294e-10;

double p7 = 1.62543e-12;

double p8 = -4.25605e-15;

double p9 = 4.67382e-18;

double y=p0+p1*pow(x,1)+p2*pow(x,2)+p3*pow(x,3)+p4*pow(x,4)+p5*pow(x,5)+p6*pow(x,6)+p7*pow(x,7)+p8*pow(x,8)+p9*pow(x,9);

return y;

}

void make_Efficiency(){

string sampleName="WIMP";

string pngFile;

string hname;

int ix=0;int iy=0; int ipad =0;

int nx=1;int ny=1;

TText *t;

gStyle->SetOptStat(kTRUE);

gStyle->SetOptStat(11111111);

gStyle->SetOptTitle(kTRUE);

f = new TFile("../betterData/H3_DataRQ_plots.root");f->cd();f->ls();

TH2D* S1Area_phd_logS2Area_phd_ss = (TH2D*) f->Get("S1Area_phd_logS2Area_phd_ss");

int nxbins = S1Area_phd_logS2Area_phd_ss->GetNbinsX();

int nybins = S1Area_phd_logS2Area_phd_ss->GetNbinsY();

int sumYAll[nxbins];

int sumYBelow[nxbins];

int sumYAbove[nxbins];

std:: cout << " opened file " << nxbins << " " << nybins << endl;

nxbins=200;

double S1bins[nxbins];

double ex[nxbins];

double Acceptance[nxbins];

double AccError[nxbins];

double Theshold[nxbins];

for(int ix=1;ix<nxbins;ix++){

int ix0=ix-1;

sumYAll[ix0]=0.;

sumYBelow[ix0]=0.;

sumYAbove[ix0]=0.;

double x = S1Area_phd_logS2Area_phd_ss->GetXaxis()->GetBinCenter(ix); //->GetBinCenterX(ix);

double MeanY = 1.0*GetMedianWimp(x);

S1bins[ix0]=x;

ex[ix]=0.0;

Theshold[ix0]=MeanY;

cout << " ix " << ix << " x " << x << " meanY " << MeanY << endl;

for(int iy=1;iy<nybins;iy++){

double y = double(S1Area_phd_logS2Area_phd_ss->GetYaxis()->GetBinCenter(iy));

double nevents = S1Area_phd_logS2Area_phd_ss->GetBinContent(ix,iy);

// cout << " ix " << ix << " iy " << iy << " y " << y << " events " << nevents << endl;

sumYAll[ix0]+=nevents;

if(y<=MeanY) {

sumYBelow[ix0]+=nevents;

}

else {

sumYAbove[ix0]+=nevents;

}

}

}

for(int ix=0;ix<nxbins-1;ix++){

cout << " ix " << ix;

cout << " all " << sumYAll[ix];

cout << " below " << sumYBelow[ix];

cout << " above " << sumYAbove[ix];

if(sumYAll[ix]>0){

double ratio = double(sumYBelow[ix])/double(sumYAll[ix]);

double error = sqrt(ratio*(1-ratio)/sumYAll[ix]);

cout << " ratio " << ratio << " +/- " << error;

Acceptance[ix]=ratio;

AccError[ix]=error;

}

cout << endl;

}

TCanvas* C1 = new TCanvas("C1","C1",1000,800);

nx=1;ny=2; std::vector<TPad*> C1Pads = make_pads(C1,nx,ny);

gStyle->SetOptStat(kFALSE);

H3Acceptance = new TGraphErrors(nxbins-1,S1bins,Acceptance,ex,AccError);

H3Acceptance = new TGraphErrors(nxbins-1,S1bins,Acceptance,ex,AccError);

GrThreshold = new TGraph(nxbins-1,S1bins,Theshold);

H3Acceptance->SetTitle("Fraction of H3 Events below the Mean of WIMP Response");

H3Acceptance->SetMarkerColor(2);

H3Acceptance->SetMarkerStyle(20);

H3Acceptance->SetMarkerSize(1);

H3Acceptance->SetMaximum(0.2);

H3Acceptance->SetMinimum(0.0);

GrThreshold->SetLineColor(2);

GrThreshold->SetLineWidth(2);

S1Area_phd_logS2Area_phd_ss->GetXaxis()->SetRangeUser(0,100);

S1Area_phd_logS2Area_phd_ss->GetYaxis()->SetRangeUser(2.0,6.0);

ix=0;iy=1; ipad = ny*ix+iy; C1->cd();C1Pads[ipad]->Draw(); C1Pads[ipad]->cd(); S1Area_phd_logS2Area_phd_ss->Draw("colz"); GrThreshold->Draw("L");

TH2D* S1AreaTemp = (TH2D*) S1Area_phd_logS2Area_phd_ss->Clone();

S1AreaTemp->Reset();

S1AreaTemp->GetXaxis()->SetRangeUser(0,100);

S1AreaTemp->GetYaxis()->SetRangeUser(0.0,0.2);

S1AreaTemp->SetTitle("Fraction of H3 Events below the Mean of WIMP Response");

ix=0;iy=0; ipad = ny*ix+iy; C1->cd();C1Pads[ipad]->Draw(); C1Pads[ipad]->cd(); S1AreaTemp->Draw(); H3Acceptance->Draw("LP");

t = new TText(40.0,0.3,"Structure in the distribution is under investigation."); t->SetTextColor(kRed+2); t->Draw();

pngFile = "H3Acceptance.png";

C1->SaveAs(pngFile.c_str());

}

Should produce:

With the data I was given:

In other words, I have had 1 meal and slept 8 hours over the past 3 days because I was sent bad data... group work is great.

void ratio_plot2(){

string hname;

int ix=0;int iy=0; int ipad =0;

int nx=1;int ny=1;

TText *t;

//get wimpl

f = new TFile("H3_DataRQ_plots.root");f->cd();f->ls();//exists

g = new TFile("WIMP_DataRQ_plots.root");g->cd();g->ls();//doesntexist

//get means of graphs and style them

TGraph *fFileMean = getMean(f, "RecoS1Area_phd_S1Bin", "RecoS2Area_phd_S1Bin");

TGraph *gFileMean = getMean(g, "RecoS1Area_phd_S1Bin", "RecoS2Area_phd_S1Bin");

fFileMean->SetLineColor(7);

gFileMean->SetLineColor(6);

fFileMean->SetLineWidth(4);

gFileMean->SetLineWidth(4);

cout<<gFileMean->GetXaxis()->GetXmax()<<endl;

//get graph from f graph

const char* graph= "S1Area_phd_logS2Area_phd_ss";

const char* reee = "wS1Area_phd_logS2Area_phd_ss";

TH2D* S1Area_phd_logS2AreaOverS1Area_phd_ss =(TH2D*) f->Get(graph); //exists

TH2D* GGS1Area_phd_logS2AreaOverS1Area_phd_ss = (TH2D*) g->Get(reee); //doesnt exist

TCanvas* C25 = new TCanvas("C25","C25",1000,1000); C25->cd();

//set axises

//S1Area_phd_logS2AreaOverS1Area_phd_ss->GetXaxis()->SetRangeUser(0,200);

//S1Area_phd_logS2AreaOverS1Area_phd_ss->GetYaxis()->SetRangeUser(0,5);

//TCanvas C_Rat* = new TCanvas("C_Rat","CRat",1000,1000); C_Rat->cd();

/*

Double_t ratioplotArr[npoints+1];

for(int i =0;i<npoints+1;i++){ratioplotArr[i]=i;}

TH1D* ratioPlot = new TH1D("ratioplot", "ratioplot ", npoints, ratioplotArr);

*/

/*not wimp mean is higher forwhatever that means

check if graph is higher than that?*/

Int_t xbins = S1Area_phd_logS2AreaOverS1Area_phd_ss->GetNbinsX();

xbins = 150; //this is the extent of thje median graph

Int_t ybins = S1Area_phd_logS2AreaOverS1Area_phd_ss->GetNbinsY();

int count = 0;

int total[xbins];

int passed[xbins];

double ratio[xbins];

for(int x =0;x<xbins;x++){

total[x]=0;passed[x]=0;ratio[x]=0.0;

}

/*

0 Bin# have it

1 BinCenter have it

2 medianWimpVal have it

3 TotalEventsInH3Sample have it

4 nonpassedevents can inverse it

5 ratio (4, 3)

6 5 error

*/

std::ofstream outputCSV;

outputCSV.open("Table.csv");

outputCSV<<"xbin, bincenter, median, totalinXbin, passed, ratio, error \n";

for(Int_t x = 0; x < xbins; x++){

Double_t delimeter = fFileMean->Eval(x);

if (x==221){

cout<<delimeter<<endl;

}

for(Int_t y = 0; y<ybins;y++){

Int_t freq = GGS1Area_phd_logS2AreaOverS1Area_phd_ss->GetBinContent(x,y);

int freqAsDouble = (int)freq;

Int_t bin = GGS1Area_phd_logS2AreaOverS1Area_phd_ss->GetBin(x,y);

Double_t someval = GGS1Area_phd_logS2AreaOverS1Area_phd_ss->GetYaxis()->GetBinCenter(y);

//cout<<someval<< " median: "<< delimeter<<endl;

if (someval>delimeter){

passed[x]+=freqAsDouble;

}

total[x]+=freqAsDouble;

}

double d = (double)delimeter;

outputCSV<<x<<","<<0<<","<<d<<","<<total[x]<<passed[x];

double ratVal = (total[x] != 0) ? passed[x]/total[x]: 0;

double err = (total[x]!=0) ? sqrt(ratVal*(1-ratVal))/total[x]: 0;

outputCSV<<","<<ratVal<<","<<err<<", \n";

}

outputCSV.close();

for(int x = 0;x<xbins;x++){

if (total[x]==0){

ratio[x] = 0;

}else{

//cout<<passed[x]<<" "<<total[x]<<endl;

ratio[x] = (double)(passed[x])/total[x];

//cout<<"x: "<<x <<": "<<ratio[x]<<endl;

}

}

S1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("colz");

//S1Area_phd_logS2AreaOverS1Area_phd_ss->SetBinContent(50,15, 100000);

GGS1Area_phd_logS2AreaOverS1Area_phd_ss->Draw("box same");

//gFileMean->Draw("same");

fFileMean->Draw("same");

C25->SaveAs("OverlayWithLines.png");

TCanvas *C_rat = new TCanvas("name", "title", 1000,1000);

C_rat->SetLogy();

Double_t ratioplotArr[xbins+1];

for(int i =0;i<xbins+1;i++){ratioplotArr[i]=i;}

TH1D* ratioPlot = new TH1D("ratioplot", "ratioplot ", xbins, ratioplotArr);

for(int i =0;i<xbins;i++){

//cout<<ratio[i]<<endl;

ratioPlot->SetBinContent(i, ratio[i]);

if (ratio[i] != 0)

ratioPlot->SetBinError(i, sqrt(ratio[i]*(1-ratio[i]))/total[i]); //meaningless error value/10

}

ratioPlot->Draw();

ratioPlot->GetYaxis()->SetRangeUser(0.001,0.1);

}

OLD STUFF!!

Log 1

cmsRel CMSSW_5_3_32

cd C(tab)

//activate virtualenv in which you can use and call cms functions

cmsenv

//make directory

mkdir Demo

//move to that directory

cd Demo

//makes a skeleton of a package containing an EDAnalyzer

    // EDAnalyzer: takes input from the event and processes the input without writing information back to the event

mkedanlzr DemoAnalyzer

//change to that directory

cd DemoAnalyzer

//builds non mutable file, skips reading the cache

scram b

cd ../..

//Run the script

cmsRun Demo/DemoAnalyzer/demoanalyzer_cfg.py

// move to source, activate the environment, call a file and view it

cd CMSSW_5_3_32/src/

cmsenv

root root://eospublic.cern.ch//eos/opendata/cms/Run2011A/ElectronHad/AOD/12Oct2013-v1/20001/001F9231-F141-E311-8F76-003048F00942.root

TBrowser t

curl http://hepcms-hn.umd.edu/~jabeen/Analysis/root-tree-example.tgz | tar xv

cd CMSSW../src/root-tree-Example

mkdir test

cd test

root -l output_VBF_HToZZTo4L_M-200_8TeV-powheg-pythia6.root

Calls structure, calls method in structure to make HiggsAnalysis class and header

HZZ4LeptonsAnalysisReduced->MakeClass("HiggsAnalysis")

Log 2

What are the fundamental particles and forces? How many particles standard model of particle physics has?

There are 4 fundamental forces, gravity, strong, weak and electromagnetism. There are six quarks, six leptons, four force-carrier particles, and the Higgs boson

What is spin of a fundamental particle? 

Spin is the degrees of rotation in which particles can exists - think back to chair demo.

What is the difference between bosons and fermions?

Fermions have half-spins, bosons have integer spins.

What is Linux?

Linux is the kernel for many different OS distributions

What is a Linux Shell?

The linux shell is a place where you can execute linux bash code

How do you find which shell your are working in?

ps -p $$

What is basic linux command structure?

Directories

What do the following commands do?

pwd, ls, ls -l, ls -lr, ls .., ls -CFx cd, cd .., mkdir, find, rmdir, which, whoami, echo, cat, more

current working directory, list items in directory and all variants, change directory, make directory, find direc/file, remove directory, what user you are, say whats in a file

echo what comes next, 

What is the difference between ls, ls ., and ls .., ls /, ls ./

Which directory your listing.

How can get help to learn more about linux commands on the command line?

google/stack overflow

What is symbol ~ used for in the directory structure?

Basically a search feature

What do uparrow and enter keys do?

Previous and future commands used.

What kind of programs/apps you can use to open the following file types: .txt, .C, .cpp, .root, .sh, .csh

Any text editor

What do operators, > and >> do?

Write to files, overwrite versus append

What is emacs?

A text editor

What is the difference between the commands

emacs -nw mycode.C 

emacs mycode.C

emacs mycode.C &

The ampersame keeps the terminal for other commands, nw means no window

What do you expect to heppen if you run the command ‘emacs myfile.txt &’ where the file myfile.txt does not exist already. 

It will make a new file

How do you save a file in emacs?

Hit the save button or  C-x C-c 

How do you search some string in emacs?

'C-s' 

LOG 2 IS IN PICTURES: I wanted all the output so I took pictures of all the output

Log 3

What is a shell script?

A program designed to run in the unix shell.

What do line #!/bin/tcsh on top of a shell script means?

Its a doctype -tells the computer the type of script to run

What is command 'chmod +x' used for?

Change permission to executable.

What do ’touch’ command do?

Makes a file

How do you execute a shell scripts?

chmod +x testscript.tcsh

./testscript

What is the argument to a shell script?

set ARG = $1

What is the difference between ’sometext’ and ‘$sometext’ in a shell script?

The dollar sign has it refer to a variable

What is piping (|) in linux and how to use it?

Chaining commands

How can you search and replace something in a file from command line or from within a shell script.

            sed -i -e 's/<find word>/<replace word>/g' test.txt

#!/bin/tcsh

#keyword ARG takes in input

set ARG=$1

# here as a bug statement we echo out the value to the screen                                                      

echo $ARG

# here we have it look for all the files in / and its subdirectories that

# contain somewhere in the file name the string we inputted

# you can change how wide or narrow your search is by changing / with ./ or ../ or /usr/ 

#finds file with argument input name in usr directory - dodes not extend                                                                                             

find /usr  -name "*$ARG*" -print

run bash file testscrippt where an argument replaces <argument>

chmod changes it to an exceutable

/testscript.tcsh <argument>

ls –l testscript.tcsh

chmod +x testscript.tcsh

ls –l testscript.tcsh

./testscript.tcsh  <argument>

makes a test file

logs homesearch output to testfile

write output of testfile

touch testfile

./homesearch.tcsh testfile

cat homesearch_testfile.txt

replaces the word one with the word of two

            sed -i -e 's/one/two/g' test.txt

log 4

Show that for small beta relativistic kinetic energy has the same expression as the usual KE (= 1/2 mv2)

Hw problem

Show that pseudorapidity and rapidity are same for massless particles.

Hw problem

Why muons travel farthest in the detector.

Hw problem

How can you create a file with emacs?

I downloaded sublime, but just call emacs in terminal

How to print a statement in c++ code.

cout<<"string"<<endl;

what is ‘ #include xyz’ mean in a C++ file? 

call some library named xyz

How to add comments in C++?

// and /**/

How do you compile c++ code to create and executable that you can then run as './executible-name’?

g++ test.cpp -o executable name

What does the g++ compiler do?

Convert c++ to machine code

What punctuation c++ code uses to separate the code lines?

';'

What does command 'g++ dumpspecs’ do?

Tells you all specs of the g++ compiler your using

How you define variables of kind int, double, float, bool, 

int x;

bool y;

float z = 1f;

double f = 10;

What do the operators ++ and —, !=, ==, do?

Add 1, subtract 1, compare not equal, compare equal

How while and for loops work in c++?

while (condition) {}

for(int i=0;i<10;i++){}

How to add debug statements in the C++code?

int debug  = 0 //debug flag

if (debug) cout<<"some statement"<<endl;

instantiates program

and prints hello world, instantiates a namespace

std::cout becomes cout

#include <iostream>

using namespace std;

int main() {

    cout <<"Hello World!" << endl;    //Print hello world to screen followed by end line (endl)

    return 0;    //Exit the program

}

prints Hello World!

compiles and runs program

• • do: g++ main.cpp

  • do: ./a.out

find /usr/include/c++ -name "4*"

/usr/include/c++/4.4.4

/usr/include/c++/4.4.7

• • g++ -dumpspecs

    • specifies the specs of the g++ compiler you have: sample output:

    • *asm:

    • %{v:-V} %{Qy:} %{!Qn:-Qy} %{n} %{T} %{Ym,*} %{Yd,*}  %{Wa,*:%*} %{m32:--32} %{m64:--64}  %{!mno-sse2avx:%{mavx:-msse2avx}} %{msse2avx:%{!mavx:-msse2avx}}

    • *asm_debug:

    • %{!g0:%{gstabs*:--gstabs}%{!gstabs*:%{g*:--gdwarf2}}} %{fdebug-prefix-map=*:--debug-prefix-map %*}

    • *asm_final:

• ls /usr/include/c++/4.4.4/

  • • algorithm           cstdarg              functional        sstream

    • array               cstdatomic           initializer_list  stack

    • backward            cstdbool             iomanip           stdatomic.h

    • bits                cstddef              ios               stdexcept

    • bitset              cstdint              iosfwd            streambuf

    • c++0x_warning.h     cstdio               iostream          string

    • cassert             cstdlib              istream           system_error

    • ccomplex            cstring              iterator          tgmath.h

    • cctype              ctgmath              limits            thread

    • cerrno              ctime                list              tr1

    • cfenv               cwchar               locale            tr1_impl

    • cfloat              cwctype              map               tuple

    • chrono              cxxabi-forced.h      memory            type_traits

    • cinttypes           cxxabi.h             mutex             typeinfo

    • ciso646             debug                new               unordered_map

    • climits             deque                numeric           unordered_set

    • clocale             exception            ostream           utility

    • cmath               exception_defines.h  parallel          valarray

    • complex             exception_ptr.h      queue             vector

    • complex.h           ext                  random            x86_64-redhat-linux

    • condition_variable  fenv.h               ratio

    • csetjmp             forward_list         regex

    • csignal             fstream              set

math.h as the only search item

• • ls /usr/include/math.h

  • /usr/include/math.h

math as a wild card tag

• • ls /usr/include/bits/*math*.h

  • /usr/include/bits/mathcalls.h   /usr/include/bits/mathinline.h

g++ --help

Options:

  -pass-exit-codes         Exit with highest error code from a phase

  --help                   Display this information

  --target-help            Display target specific command line options

  --help={target|optimizers|warnings|params|[^]{joined|separate|undocumented}}[,...]

                           Display specific types of command line options

  (Use '-v --help' to display command line options of sub-processes)

  --version                Display compiler version information

  -dumpspecs               Display all of the built in spec strings

  -dumpversion             Display the version of the compiler

  -dumpmachine             Display the compiler's target processor

  -print-search-dirs       Display the directories in the compiler's search path

  -print-libgcc-file-name  Display the name of the compiler's companion library

  -print-file-name=<lib>   Display the full path to library <lib>

  -print-prog-name=<prog>  Display the full path to compiler component <prog>

  -print-multi-directory   Display the root directory for versions of libgcc

  -print-multi-lib         Display the mapping between command line options and

                           multiple library search directories

save and choose executable name

g++ -O mypreferredname.out mycode.C

#include <iostream>

using namespace std;

//all executable code needs to be called in main

int main() {

        //print hello world

    cout << "hello world" << end;

    int i=2;

    //print i

    cout << "i = " <<i<<endl;

    double a=3.3;

    //double type

    cout << "a = " <<a<<endl;

    int j = a*i;

    cout << "a*i = "<<j<<endl;

    return 0;

}

out: hello world i = 2 a = 3.3 a*i = 6

#include <iostream>

using namespace std;

int main() {

    int n=10;

    cout << "n is "<<n<<endl;

     //decrement value

    n--;

    cout<<"n is now "<<n<<endl;

     //increment value

    n++;

    cout<<n is now "<<n<<endl;

    return 0;

}

out: . n is 10 n is now 9 n is now 10

++n versus n++, ignoring efficiency, the difference is that in ++n the ++ has precedence before the variable is used in any other in line expression while the appended operand version does not have precednece

#include <iostream>

using namespace std;

int main() {

    bool prop;

        //ture statement

    prop = (5>1);

    cout<<"prop is "<<prop<<endl;

     //false statement

    prop = (1>5);

    cout<<"prop is "<<prop<<endl;

     //true statement

    prop = (1 != 5);

    cout << "prop is " <<prop<<endl;

    return 0;

}

out: prop is 1 prop is 0 prop is 1

#include <iostream>

using namespace std;

int main() {

    int n=10;

       //while loop, decremention loop

    while(n>0) {

        cout<<"n is "<<n<<endl;

        n--;

    }

    return 0;

}

out: n is 10 n is 9 n is 8 n is 7 n is 6 n is 5 n is 4 n is 3 n is 2 n is 1

#include <iostream>

using namespace std;

int main() {

    // when we declare a for loop, we also initialize the loop variable,

    // specify the exit condition, and tell the program how to modify the

    // loop variable at the end of each loop

    for (int n=10; n>0; n--) {

        cout<<"n is "<<n<<endl;

    }

    // in a for loop, the loop variable (in this case, 'n') only exists in

    // the loop. we are not able to call 'n' from out here

    // uncomment the following line and see for yourself

    // cout<<"n outside the loop: "<<n;

    return 0;

}

out: n is 10 n is 9 n is 8 n is 7 n is 6 n is 5 n is 4 n is 3 n is 2 n is 1

#include <iostream>

using namespace std;

int main() {

    int n=0, m=0;

      // double while loop - unrollable why would you do htis

    while(n<10) {

        // this is the slow (or outer) loop

        cout << "n is " << n << ": ";

        m=0;

           while(m<=n) {

            // this is the fast (or inner) loop

            // in this loop, the slow loop variable (n) is a constant

            // this loop must run to completion before the slow loop

            // can progress (during every iteration of the slow loop!)

            cout << m;

            m++;

        }

        // now the fast loop has finished and the slow loop can

        // continue with the current iteration

        cout << endl;

        n++;

    }

    return 0 ;

}  

out: n is 0: 0 n is 1: 01 n is 2: 012 n is 3: 0123 n is 4: 01234 n is 5: 012345 n is 6: 0123456 n is 7: 01234567 n is 8: 012345678 n is 9: 0123456789

Previous program with debug:

the current value of my variable 5 n is 10 n is 9 n is 8 n is 7 n is 6 n is 5 n is 4 n is 3 n is 2 n is 1 

Ideas of debgug:

flag at top of program: i.e int Debug = 0;

then on check statemements: if (Debug) cout<<"Message"<<endl;

GDB debugging:

g++ -g file.cpp //compile with debug flag

gdb //run gdb

file a.out // call executable in gdb

list //list some of the files contents

list 12 //list some file contents centered at 12

break 12 //make a break at line 12

run //run until break

info locals //get info of local variables

c //continue program normally

step //step through program step by step

It is important to note GDB debugs in set states which can be useful, George Hotz made Qira which presents a kind of stateless debugging, This can be useful as you look for all instances a pointer was effected

which can make for more useful debugging.

https://qira.me/

#include <iostream>

using namespace std;

int main() {

  int i = 1;

  while (i<1.5e9) {

    i *= 10;

  }

  cout << "The final value of 'i' was: " << i << endl;

  return 0;

}

This will break, before it prints anything: as you'll have an integer overflow error, both a long or a double could hold the value size necessary.

log 5

How to read a grapgh and extract values for the hysics quantities involved?

Calculate tranverse momentum of a particle moving in electric and magnetic field.

Logic statements in C++,  if-else.

if (condition) {}else{}

How you define Pointers in C++? How you can get the address and the value of a pointer?

int *p = new int(5);

cout<<p<<endl;

if non pointer: &i;

How can you define a pointer with a new memory address?

int *p = new int(5);

Using ‘date' and ‘awk' commands in a script.

#include <ctime>

time_t now = time(0)

char* dt = ctime(&now);

#include <iostream>

using namespace std;

int main() {

     int i = 10

       // reference memory address of i with ampersame

     cout << "The memory address of i is " << &i << "\n";

     cout << "The data stored at memory address " << &i << " is " << i << "\n";

        //set pointer p to the memory address of I

     int* p = &i;

     cout << "The value of p is " << p << "\n";

     cout << "We say that p 'points at' the memory location referenced by address " << p << "\n";

      //Dereference pointer with asterisk

     cout << "The data stored at memory address " << p << " is " << *p << "\n";

     return 0;

}

prints the value of p, the pointer memory address, repeatedly..

/*PROGRAM 1*/

#include <iostream>

using namespace std;

int main(){

    //instantiate two variables, same value, different memory address

     int i = 10;

     int j = i;

     cout << "i= " << i << " and j= " << j << "\n";

     i=5;

     cout << "i= " << i << " and j= " << j << "\n";

     j=1;

     cout << "i= " << i << " and j= " << j << "\n";

     return 0;

}

final variable states:

i =5

j =1

/*PROGRAM 2*/

#include <iostream>

using namespace std;

int main(){

    //instantiate normal integer type primitive variable

     int i = 10;

        //make second variable p, a pointer which references the integer type and set the pointer to the same location in memory as i

     int* p = &i;

     cout << "i= " << i << " and *p= " << *p << "\n";

     i=5;

     cout << "i= " << i << " and *p= " << *p << "\n";

     *p=1;

     cout << "i= " << i << " and *p= " << *p << "\n";

     return 0;

}

final variable states:

i = 1

p = 1

While the end behavior of both these programs are similar, the second program is more efficient. The first program instantiates 2 distinct pointers which both contain the value 10. The second program instantiates one pointer than creates a second variable which references the same pointer which is more memory efficient. In addition, changing one variable in the second program will effect both, while the variables are completely independent in program 1.

#include <iostream>

using namespace std;

int main(){

     int* p = new int(5);

     cout << "p points at address " << p << "\n";

     cout << "The data stored in address " << p << " is " << *p << "\n";

     *p = 10;

     cout << "Now the data stored in address " << p << " is " << *p << "\n";

     return 0;

}

Here, the new operator allows the program to instantiate a pointer type of int without having to reference an already created variable. "new" saves the 5 in a unique memory address and allows the programmer to reference it. 

#!/bin/tcsh

//gets current date

set month = `date | awk '{print $2}'`

    echo $month

set tmpday = `date | awk '{print $3}'`

    echo $tmpday

//why is this here, why are you correcting date

if ($tmpday < 10) then

    set today =  " $tmpday"

else

    set today = $tmpday

endif

echo $today

//prints the files

ls -lat | grep "$month $today"

prints date: then all files modified on that date and info on said file

Oct

8

8

drwxrwxr-x 2 cms-opendata cms-opendata  4096 Oct  8 12:10 .

-rwxrwxr-x 1 cms-opendata cms-opendata   255 Oct  8 12:10 l5.tcsh

-rwxrwxr-x 1 cms-opendata cms-opendata   254 Oct  8 12:08 log5.tcsh

-rw-rw-r-- 1 cms-opendata cms-opendata    13 Oct  8 12:03 vectors.txt

-rw-rw-r-- 1 cms-opendata cms-opendata  1221 Oct  8 12:01 useDotprod.cpp

-rwxrwxr-x 1 cms-opendata cms-opendata 10576 Oct  8 11:56 a.out

-rw-rw-r-- 1 cms-opendata cms-opendata   213 Oct  8 11:55 dotprod.cpp

-rw-rw-r-- 1 cms-opendata cms-opendata    16 Oct  8 11:52 example.txt

-rw-rw-r-- 1 cms-opendata cms-opendata   200 Oct  8 11:51 read1.cpp

-rw-rw-r-- 1 cms-opendata cms-opendata   543 Oct  8 11:49 arr2.cpp

-rw-rw-r-- 1 cms-opendata cms-opendata   425 Oct  8 11:44 arr1.cpp

drwxr-xr-x 7 cms-opendata cms-opendata  4096 Oct  8 11:43 ..

Downloading sublime:

wget https://download.sublimetext.com/sublime_text_3_build_3126_x64.tar.bz2

sudo tar -vxjf sublime_text_3_build_3126_x64.tar.bz2 -C /opt

sudo ln -s /opt/sublime_text_3/sublime_text /usr/bin/sublime3

Log 6

#include <iostream>

using namespace std;

int main() {

        //array instantiation

    int ii[3] = {1,2,3};

    int j=0;

    while (j<3) {

        cout <<" ii of "<<j<<" is "<<ii[j]<<endl;

        j++;

    }

     //horrible array instantiation

    int LL[2][3] = {1,2,3,4,5,6};

    j=0;

    int k;

    while (j<2) {

        k=0;

        while (k<3) {

            cout<<" LL of "<<j<<" "<<k<<" is "<<LL[j][k]<<endl;

            k++;

        }

        j++;

    }

    return 0;

}

ii of 0 is 1

 ii of 1 is 2

 ii of 2 is 3

 LL of 0 0 is 1

 LL of 0 1 is 2

 LL of 0 2 is 3

 LL of 1 0 is 4

 LL of 1 1 is 5

 LL of 1 2 is 6

//literal copy of above

#include <iostream>

using namespace std;

int main() {

    // a loop to demonstrate 1D arrays

    int ii[3] = {1,2,3};

    int j=0;

    while (j<3) {

        cout <<" ii of "<<j<<" is "<<ii[j]<<endl;

        j++;

    }

    // a loop ot demonstrate 2D arrays

    int LL[2][3] = {1,2,3,4,5,6};

    j=0;

    int k;

    while (j<2) {

        k=0;    // do not forget to initialize k here

        while (k<3) {

            cout<<" LL of "<<j<<" "<<k<<" is "<<LL[j][k]<<endl;

            k++;

        }

        j++;

    }

    return 0;

}

ii of 0 is 1

 ii of 1 is 2

 ii of 2 is 3

 LL of 0 0 is 1

 LL of 0 1 is 2

 LL of 0 2 is 3

 LL of 1 0 is 4

 LL of 1 1 is 5

 LL of 1 2 is 6

#include <iostream>

#include <fstream>

using namespace std;

int main() {

    ofstream myfile;

        //opens sample text file

    myfile.open("example.txt");

     //writes to text file

    myfile<<"write some junk.";

        //memory leaks are good

    myfile.close();

    return 0;

}

Writes "write some junk" to the file example.txt

CALLED DOTPROD FILE

#include <iostream>

using namespace std;

//double function that only works with 3D vectors

double dot_prod(double v1[3],double v2[3]) {

  double dotdot;

    //element wise multiplication

  dotdot = v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];

  cout<<" The dot product is "<<dotdot<<endl;

  return 0;

}

//use dot prod file

//reads from input file, does dot product

#include <iostream>

#include <fstream>

// include the program dotprod.cpp so that we can find the dot_prod function                                         

#include "dotprod.cpp"

using namespace std;

int main () {

  // declare the vectors                                                                                             

  double vector1[3];

  double vector2[3];

  // open the input file                                                                                             

  ifstream infile;

  infile.open("vectors.txt");

//sets vector1 to the firxt line vector

  // store the input in the vectors and print the vectors for the user                                               

  infile>>vector1[0]>>vector1[1]>>vector1[2];

  cout<<" Vector 1 is ("<<vector1[0]<<","<<vector1[1]<<","<<vector1[2]<<")"<<endl;

//sets vec2 to second line

  infile>>vector2[0]>>vector2[1]>>vector2[2];

  cout<<" Vector 2 is ("<<vector2[0]<<","<<vector2[1]<<","<<vector2[2]<<")"<<endl;

  // close the input file                                                                                            

  infile.close();

  // call the dot_prod function from dotprod.cpp  

    //calls dotprod function from dotprod file                                                                   

  dot_prod(vector1,vector2);

  return 0;

}

Vector 1 is (1,2,3)

 Vector 2 is (4,5,6)

 The dot product is 32

Code for scalar dotProd.

#include <iostream>

using namespace std;

//double function that only works with 3D vectors

double scaleProd(double scalar,double v2[3]) {

  double dotdot;

    //element wise multiplication

  v[0]*=scalar;

    v[1]*=scalar;

    v[2]*=scalar

  cout<<" Vector 2 is ("<<v2[0]<<","<<v2[1]<<","<<v2[2]<<")"<<endl;

  return 0;

}

#include <iostream>

#include <math.h>

using namespace std;

// this function is the actual random number generator

// this code is stolen from the book numerical recipes for fortran

// it relies on random generated from overflow of memory locations

// and is a pseudo random number generator

const int a = 7141;

const int c = 54773;

const int mmod=256200;

//reads a value as a pointer thats pretty damn clever

//magical operations to get new random number

double getFlatRandom(int& inew) {

double mranflat = 0.;

//mod it for overflow reasons

inew = inew%mmod;

double aa = double(inew)/double(mmod);

mranflat=aa;

//make it larger for reasons?

inew = a*inew+c;

return mranflat;

}

// in this code, we will call the pseudo-random number generator and learn some things

// about its properties by filling and then displaying a histogram

int main() {

    int num;

    cout << "Enter the number of loop iterations: ";

    cin >> num;

int inew = 2345;  // This is the "seed" for the random number generator

// we will put the results from the call into a histogram that we can look at, to learn some of its

// properties. This histogram has 10 bins. 

int histo[10] = {0,0,0,0,0,0,0,0,0,0};

double atmp; // 

// call the random number generator 1000 times and fill a histogram

for (int i = 0; i<num; i++){

                //passing in an integer value as a pointer is pretty damn clever

         atmp=getFlatRandom(inew);  // call the random number generator

histo[int(atmp*10)]++;  // increment the histogram bin the number falls within

}

// print the histogram to the screen

for (int i = 0; i<10; i++){

cout<<i<<": ";

for (int j=0; j<int((double(100*histo[i])/double(num))+0.5); j++) {

            cout << "=";

        }

    cout << endl;

        }

return 0;

}

Log 7

Edited code - boolean made to boolean, N to 1000 and added closing } to for loop

#include "Riostream.h"

#include <iostream>

#include <fstream>

/* 

some code to help us understand why and how resolutions affect our ability to measure the

mass of a particle

*/

// main routine

void resolutions(){

// set debug mode

bool idebug=false;

// set up the display

gROOT->Reset();

gROOT->SetStyle("Plain");

gStyle->SetOptStat(0);

gStyle->SetOptFit(11);

gStyle->SetErrorX(0);

c1 = new TCanvas("c1","Our data",200,10,700,500);

c1->SetFillColor(0);

c1->SetBorderMode(0);

c1->SetBorderSize(1);

c1->SetFrameBorderMode(0);

// book some histograms

TH1F *histo1 = new TH1F("histo1","mass",1000,0.,200.);

// setup random number generator

gRandom->SetSeed();

// get the secret parameters

ifstream myfile;

myfile.open("secretparameters.txt");

double mass, resE, resA;

myfile>>mass>>resE>>resA;

myfile.close();

// make N bosons at rest

// each boson will decay to 2 back-to-back daughter particles

// they will be back-to-back by conservation of momentum, since the momentum of the mother was zero

// assume that the mass of the daughter particles is small compared to the mother mass so we can

// assume that their energy will be large compared to their mass and we can treat them as massless.

// and thus their energy is equal to the magnitude of their momenta

// by conservations of energy, their energy

// work in 2D. Can always choose my coordinate system so that the 2 daughters are in the same plane

int N=1000;

double etrue1,etrue2,phitrue1,phitrue2; // true energy of the 2 daughters

double e1,px1,py1,phi1; // smeared 4-momenta of daughter 1

double e2,px2,py2,phi2; // smeared 4-momenta of daughter 2

double masssmeared;

for(int i=0;i<N;i++) {

// set true energy

etrue1=mass/2.;

etrue2=mass/2.;

if(idebug) cout<<"etrue "<<etrue1<<" "<<etrue2<<endl;

// choose phi for daughter 1 and daughter 2

phitrue1=2*TMath::Pi()*gRandom->Rndm();

phitrue2 = phitrue1+TMath::Pi();

if(idebug) cout<<"phitrue "<<phitrue1<<" "<<phitrue2<<endl;

// smear true energy with resolution of detector

e1=etrue1+resE*gRandom->Gaus(0.,1.);

e2=etrue2+resE*gRandom->Gaus(0.,1.);

if(idebug) cout<<"e "<<e1<<" "<<e2<<endl;

//smear angles with resolution of the detector

phi1=phitrue1+resA*gRandom->Gaus(0.,1.);

phi2=phitrue2+resA*gRandom->Gaus(0.,1.);

if(idebug) cout<<"phi "<<phi1<<" "<<phi2<<endl;

//calculate 4 momenta after smearing

px1=e1*cos(phi1);

py1=e1*sin(phi1);

px2=e2*cos(phi2);

py2=e2*sin(phi2);

if(idebug) cout<<"pxs "<<px1<<" "<<py1<<" "<<px2<<" "<<py2<<endl;

// calculate smeared mass

masssmeared=sqrt((e1+e2)*(e1+e2) - (px1+px2)*(px1+px2) - (py1+py2)*(py1+py2));

if(idebug) cout<<"masssmeared "<<masssmeared<<endl;

histo1->Fill(masssmeared);

histo1->Draw("");}

c1->Update();

c1->SaveAs("c1.gif");

}

}

full code for image lager - after madgraph

MADGRAPH

[cms-opendata@localhost CMSSW_5_3_32]$ ./bin/mg5_MC

bash: ./bin/mg5_MC: No such file or directory

[cms-opendata@localhost CMSSW_5_3_32]$ /bin/mg5_MC

bash: /bin/mg5_MC: No such file or directory

[cms-opendata@localhost CMSSW_5_3_32]$ cd MG5_aMC_v2_6_7     

[cms-opendata@localhost MG5_aMC_v2_6_7]$ ls

HELAS    MadSpin  Template         aloha  doc.tgz   mg5decay       tests

INSTALL  PLUGIN   UpdateNotes.txt  bin    input     models         vendor

LICENSE  README   VERSION          doc    madgraph  proc_card.dat

[cms-opendata@localhost MG5_aMC_v2_6_7]$ bin/mg5_aMC 

************************************************************

*                                                          *

*                     W E L C O M E to                     *

*              M A D G R A P H 5 _ a M C @ N L O           *

*                                                          *

*                                                          *

*                 *                       *                *

*                   *        * *        *                  *

*                     * * * * 5 * * * *                    *

*                   *        * *        *                  *

*                 *                       *                *

*                                                          *

*         VERSION 2.6.7                 2019-10-16         *

*                                                          *

*    The MadGraph5_aMC@NLO Development Team - Find us at   *

*    https://server06.fynu.ucl.ac.be/projects/madgraph     *

*                            and                           *

*            http://amcatnlo.web.cern.ch/amcatnlo/         *

*                                                          *

*               Type 'help' for in-line help.              *

*           Type 'tutorial' to learn how MG5 works         *

*    Type 'tutorial aMCatNLO' to learn how aMC@NLO works   *

*    Type 'tutorial MadLoop' to learn how MadLoop works    *

*                                                          *

************************************************************

load MG5 configuration from input/mg5_configuration.txt 

python2.6 does not support such functionalities please use python2.7

Option low_mem_multicore_nlo_generation from config file not understood

fastjet-config does not seem to correspond to a valid fastjet-config executable (v3+). We will use fjcore instead.

 Please set the 'fastjet'variable to the full (absolute) /PATH/TO/fastjet-config (including fastjet-config).

 MG5_aMC> set fastjet /PATH/TO/fastjet-config

lhapdf-config does not seem to correspond to a valid lhapdf-config executable. 

Please set the 'lhapdf' variable to the (absolute) /PATH/TO/lhapdf-config (including lhapdf-config).

Note that you can still compile and run aMC@NLO with the built-in PDFs

 MG5_aMC> set lhapdf /PATH/TO/lhapdf-config

Using default text editor "vi". Set another one in ./input/mg5_configuration.txt

Using default eps viewer "evince". Set another one in ./input/mg5_configuration.txt

Using default web browser "firefox". Set another one in ./input/mg5_configuration.txt

Loading default model: sm

INFO: load particles 

INFO: load vertices 

INFO: Restrict model sm with file models/sm/restrict_default.dat . 

INFO: Run "set stdout_level DEBUG" before import for more information. 

INFO: Change particles name to pass to MG5 convention 

Defined multiparticle p = g u c d s u~ c~ d~ s~

Defined multiparticle j = g u c d s u~ c~ d~ s~

Defined multiparticle l+ = e+ mu+

Defined multiparticle l- = e- mu-

Defined multiparticle vl = ve vm vt

Defined multiparticle vl~ = ve~ vm~ vt~

Defined multiparticle all = g u c d s u~ c~ d~ s~ a ve vm vt e- mu- ve~ vm~ vt~ e+ mu+ t b t~ b~ z w+ h w- ta- ta+

MG5_aMC>generate p p >t t~

Command "generate p p >t t~" interrupted with error:

InvalidCmd : No particle >t in model

MG5_aMC>generate p p > t t~

INFO: Checking for minimal orders which gives processes. 

INFO: Please specify coupling orders to bypass this step. 

INFO: Trying coupling order WEIGHTED<=2: WEIGTHED IS 2*QED+QCD 

INFO: Trying process: g g > t t~ WEIGHTED<=2 @1  

INFO: Process has 3 diagrams 

INFO: Trying process: u u~ > t t~ WEIGHTED<=2 @1  

INFO: Process has 1 diagrams 

INFO: Trying process: u c~ > t t~ WEIGHTED<=2 @1  

INFO: Trying process: c u~ > t t~ WEIGHTED<=2 @1  

INFO: Trying process: c c~ > t t~ WEIGHTED<=2 @1  

INFO: Process has 1 diagrams 

INFO: Trying process: d d~ > t t~ WEIGHTED<=2 @1  

INFO: Process has 1 diagrams 

INFO: Trying process: d s~ > t t~ WEIGHTED<=2 @1  

INFO: Trying process: s d~ > t t~ WEIGHTED<=2 @1  

INFO: Trying process: s s~ > t t~ WEIGHTED<=2 @1  

INFO: Process has 1 diagrams 

INFO: Process u~ u > t t~ added to mirror process u u~ > t t~ 

INFO: Process c~ c > t t~ added to mirror process c c~ > t t~ 

INFO: Process d~ d > t t~ added to mirror process d d~ > t t~ 

INFO: Process s~ s > t t~ added to mirror process s s~ > t t~ 

5 processes with 7 diagrams generated in 0.081 s

Total: 5 processes with 7 diagrams

MG5_aMC>display processes

Process: g g > t t~ WEIGHTED<=2 @1

Process: u u~ > t t~ WEIGHTED<=2 @1

Process: c c~ > t t~ WEIGHTED<=2 @1

Process: d d~ > t t~ WEIGHTED<=2 @1

Process: s s~ > t t~ WEIGHTED<=2 @1

MG5_aMC>display particles

Current model contains 17 particles:

w+/w- ve/ve~ vm/vm~ vt/vt~ u/u~ c/c~ t/t~ d/d~ s/s~ b/b~ e-/e+ mu-/mu+ ta-/ta+ 

a z g h 

MG5_aMC>display multiparticles

Multiparticle labels:

all = g u c d s u~ c~ d~ s~ a ve vm vt e- mu- ve~ vm~ vt~ e+ mu+ t b t~ b~ z w+ h w- ta- ta+

l- = e- mu-

j = g u c d s u~ c~ d~ s~

vl = ve vm vt

l+ = e+ mu+

p = g u c d s u~ c~ d~ s~

vl~ = ve~ vm~ vt~

MG5_aMC>

MG5_aMC>add process p p > W + j, W+ > l+ vl

Command "add process p p > W + j, W+ > l+ vl" interrupted with error:

InvalidCmd : No particle w in model

MG5_aMC>add process p p > W+j, W+ > l+ vl

Command "add process p p > W+j, W+ > l+ vl" interrupted with error:

InvalidCmd : No particle w+j in model

MG5_aMC>add process p p > W+ j, W+ > l+ vl

INFO: Checking for minimal orders which gives processes. 

INFO: Please specify coupling orders to bypass this step. 

INFO: Trying coupling order WEIGHTED<=3: WEIGTHED IS 2*QED+QCD 

INFO: Trying process: g u > w+ d WEIGHTED<=3 @2  

INFO: Process has 2 diagrams 

INFO: Trying process: g u > w+ s WEIGHTED<=3 @2  

INFO: Trying process: g c > w+ d WEIGHTED<=3 @2  

INFO: Trying process: g c > w+ s WEIGHTED<=3 @2  

INFO: Process has 2 diagrams 

INFO: Crossed process found for g d~ > w+ u~, reuse diagrams. 

INFO: Crossed process found for g s~ > w+ c~, reuse diagrams. 

INFO: Process u g > w+ d added to mirror process g u > w+ d 

INFO: Crossed process found for u d~ > w+ g, reuse diagrams. 

INFO: Process c g > w+ s added to mirror process g c > w+ s 

INFO: Crossed process found for c s~ > w+ g, reuse diagrams. 

INFO: Process d~ g > w+ u~ added to mirror process g d~ > w+ u~ 

INFO: Process d~ u > w+ g added to mirror process u d~ > w+ g 

INFO: Process s~ g > w+ c~ added to mirror process g s~ > w+ c~ 

INFO: Process s~ c > w+ g added to mirror process c s~ > w+ g 

INFO: Checking for minimal orders which gives processes. 

INFO: Please specify coupling orders to bypass this step. 

INFO: Trying process: w+ > e+ ve WEIGHTED<=2  

INFO: Process has 1 diagrams 

INFO: Trying process: w+ > e+ vm WEIGHTED<=2  

INFO: Trying process: w+ > e+ vt WEIGHTED<=2  

INFO: Trying process: w+ > mu+ ve WEIGHTED<=2  

INFO: Trying process: w+ > mu+ vm WEIGHTED<=2  

INFO: Process has 1 diagrams 

INFO: Trying process: w+ > mu+ vt WEIGHTED<=2  

1 processes with 14 diagrams generated in 0.186 s

Total: 6 processes with 21 diagrams

MG5_aMC>output MY_FIRST_MG5_RUN

INFO: initialize a new directory: MY_FIRST_MG5_RUN 

INFO: remove old information in MY_FIRST_MG5_RUN 

INFO: Organizing processes into subprocess groups 

INFO: Organizing processes into subprocess groups 

INFO: Organizing processes into subprocess groups 

INFO: Generating Helas calls for process: g u > w+ d WEIGHTED<=3 @2 

INFO: Combined process g c > w+ s WEIGHTED<=3 @2 with process g u > w+ d WEIGHTED<=3 @2 

INFO: Generating Helas calls for process: g d~ > w+ u~ WEIGHTED<=3 @2 

INFO: Combined process g s~ > w+ c~ WEIGHTED<=3 @2 with process g d~ > w+ u~ WEIGHTED<=3 @2 

INFO: Generating Helas calls for process: u d~ > w+ g WEIGHTED<=3 @2 

INFO: Combined process c s~ > w+ g WEIGHTED<=3 @2 with process u d~ > w+ g WEIGHTED<=3 @2 

INFO: Generating Helas calls for process: w+ > e+ ve WEIGHTED<=2 

INFO: Generating Helas calls for process: w+ > mu+ vm WEIGHTED<=2 

INFO: Combine g u > w+ d WEIGHTED<=3 @2 with decays w+ > e+ ve WEIGHTED<=2 

INFO: Combine g u > w+ d WEIGHTED<=3 @2 with decays w+ > mu+ vm WEIGHTED<=2 

INFO: Combining process with g u > w+ d WEIGHTED<=3 @2

  Decay: w+ > e+ ve WEIGHTED<=2 

INFO: Combine g d~ > w+ u~ WEIGHTED<=3 @2 with decays w+ > e+ ve WEIGHTED<=2 

INFO: Combine g d~ > w+ u~ WEIGHTED<=3 @2 with decays w+ > mu+ vm WEIGHTED<=2 

INFO: Combining process with g d~ > w+ u~ WEIGHTED<=3 @2

  Decay: w+ > e+ ve WEIGHTED<=2 

INFO: Combine u d~ > w+ g WEIGHTED<=3 @2 with decays w+ > e+ ve WEIGHTED<=2 

INFO: Combine u d~ > w+ g WEIGHTED<=3 @2 with decays w+ > mu+ vm WEIGHTED<=2 

INFO: Combining process with u d~ > w+ g WEIGHTED<=3 @2

  Decay: w+ > e+ ve WEIGHTED<=2 

INFO: Processing color information for process: g u > w+ d @2

  Decay: w+ > e+ ve WEIGHTED<=2 

INFO: Processing color information for process: g d~ > w+ u~ @2

  Decay: w+ > e+ ve WEIGHTED<=2 

INFO: Processing color information for process: u d~ > w+ g @2

  Decay: w+ > e+ ve WEIGHTED<=2 

INFO: Generating Helas calls for process: g g > t t~ WEIGHTED<=2 @1 

INFO: Processing color information for process: g g > t t~ @1 

INFO: Generating Helas calls for process: u u~ > t t~ WEIGHTED<=2 @1 

INFO: Processing color information for process: u u~ > t t~ @1 

INFO: Combined process c c~ > t t~ WEIGHTED<=2 @1 with process u u~ > t t~ WEIGHTED<=2 @1 

INFO: Combined process d d~ > t t~ WEIGHTED<=2 @1 with process u u~ > t t~ WEIGHTED<=2 @1 

INFO: Combined process s s~ > t t~ WEIGHTED<=2 @1 with process u u~ > t t~ WEIGHTED<=2 @1 

INFO: Creating files in directory P1_gg_ttx 

INFO: Generating Feynman diagrams for Process: g g > t t~ WEIGHTED<=2 @1 

INFO: Finding symmetric diagrams for subprocess group gg_ttx 

INFO: Creating files in directory P1_qq_ttx 

INFO: Generating Feynman diagrams for Process: u u~ > t t~ WEIGHTED<=2 @1 

INFO: Finding symmetric diagrams for subprocess group qq_ttx 

INFO: Creating files in directory P2_gq_wpq_wp_lvl 

INFO: Generating Feynman diagrams for Process: g u > w+ d WEIGHTED<=3 @2

  Decay: w+ > e+ ve WEIGHTED<=2 

INFO: Generating Feynman diagrams for Process: g d~ > w+ u~ WEIGHTED<=3 @2

  Decay: w+ > e+ ve WEIGHTED<=2 

INFO: Finding symmetric diagrams for subprocess group gq_wpq_wp_lvl 

INFO: Creating files in directory P2_qq_wpg_wp_lvl 

INFO: Generating Feynman diagrams for Process: u d~ > w+ g WEIGHTED<=3 @2

  Decay: w+ > e+ ve WEIGHTED<=2 

INFO: Finding symmetric diagrams for subprocess group qq_wpg_wp_lvl 

Generated helas calls for 5 subprocesses (10 diagrams) in 0.147 s

Wrote files for 46 helas calls in 0.372 s

ALOHA: aloha creates FFV2 routines

ALOHA: aloha creates FFV1 routines

ALOHA: aloha creates VVV1 set of routines with options: P0

save configuration file to /home/cms-opendata/CMSSW_5_3_32/MG5_aMC_v2_6_7/MY_FIRST_MG5_RUN/Cards/me5_configuration.txt

INFO: Use Fortran compiler gfortran 

INFO: Use c++ compiler g++ 

INFO: Generate jpeg diagrams 

INFO: Generate web pages 

Output to directory /home/cms-opendata/CMSSW_5_3_32/MG5_aMC_v2_6_7/MY_FIRST_MG5_RUN done.

Type "launch" to generate events from this process, or see

/home/cms-opendata/CMSSW_5_3_32/MG5_aMC_v2_6_7/MY_FIRST_MG5_RUN/README

Run "open index.html" to see more information about this process.

MG5_aMC>launch MY_FIRST_MG5_RUN/

************************************************************

*                                                          *

*                      W E L C O M E to                    *

*             M A D G R A P H 5 _ a M C @ N L O            *

*                      M A D E V E N T                     *

*                                                          *

*                 *                       *                *

*                   *        * *        *                  *

*                     * * * * 5 * * * *                    *

*                   *        * *        *                  *

*                 *                       *                *

*                                                          *

*         VERSION 2.6.7                 2019-10-16         *

*                                                          *

*    The MadGraph5_aMC@NLO Development Team - Find us at   *

*    https://server06.fynu.ucl.ac.be/projects/madgraph     *

*                                                          *

*               Type 'help' for in-line help.              *

*                                                          *

************************************************************

INFO: load configuration from /home/cms-opendata/CMSSW_5_3_32/MG5_aMC_v2_6_7/MY_FIRST_MG5_RUN/Cards/me5_configuration.txt  

INFO: load configuration from /home/cms-opendata/CMSSW_5_3_32/MG5_aMC_v2_6_7/input/mg5_configuration.txt  

INFO: load configuration from /home/cms-opendata/CMSSW_5_3_32/MG5_aMC_v2_6_7/MY_FIRST_MG5_RUN/Cards/me5_configuration.txt  

Using default text editor "vi". Set another one in ./input/mg5_configuration.txt

generate_events run_01

The following switches determine which programs are run:

/===========================================================================\

| 1. Choose the shower/hadronization program     shower = Not Avail.        |

| 2. Choose the detector simulation program    detector = Not Avail.        |

| 3. Choose an analysis package (plot/convert) analysis = Not Avail.        |

| 4. Decay onshell particles                    madspin = OFF               |

| 5. Add weights to events for new hypp.       reweight = OFF               |

\===========================================================================/

Either type the switch number (1 to 5) to change its setting,

Set any switch explicitly (e.g. type 'madspin=ON' at the prompt)

Type 'help' for the list of all valid option

Type '0', 'auto', 'done' or just press enter when you are done.[60s to answer] 

>launch MY_FIRST_MG5_RUN 

Not valid command: launch MY_FIRST_MG5_RUN

The following switches determine which programs are run:

/===========================================================================\

| 1. Choose the shower/hadronization program     shower = Not Avail.        |

| 2. Choose the detector simulation program    detector = Not Avail.        |

| 3. Choose an analysis package (plot/convert) analysis = Not Avail.        |

| 4. Decay onshell particles                    madspin = OFF               |

| 5. Add weights to events for new hypp.       reweight = OFF               |

\===========================================================================/

Either type the switch number (1 to 5) to change its setting,

Set any switch explicitly (e.g. type 'madspin=ON' at the prompt)

Type 'help' for the list of all valid option

Type '0', 'auto', 'done' or just press enter when you are done.

The following switches determine which programs are run:

/===========================================================================\

| 1. Choose the shower/hadronization program     shower = Not Avail.        |

| 2. Choose the detector simulation program    detector = Not Avail.        |

| 3. Choose an analysis package (plot/convert) analysis = Not Avail.        |

| 4. Decay onshell particles                    madspin = OFF               |

| 5. Add weights to events for new hypp.       reweight = OFF               |

\===========================================================================/

Either type the switch number (1 to 5) to change its setting,

Set any switch explicitly (e.g. type 'madspin=ON' at the prompt)

Type 'help' for the list of all valid option

Type '0', 'auto', 'done' or just press enter when you are done.

>0

Do you want to edit a card (press enter to bypass editing)?

/------------------------------------------------------------\

|  1. param : param_card.dat                                 |

|  2. run   : run_card.dat                                   |

\------------------------------------------------------------/

 you can also

   - enter the path to a valid card or banner.

   - use the 'set' command to modify a parameter directly.

     The set option works only for param_card and run_card.

     Type 'help set' for more information on this command.

   - call an external program (ASperGE/MadWidth/...).

     Type 'help' for the list of available command

 [0, done, 1, param, 2, run, enter path][90s to answer] 

>2

open /home/cms-opendata/CMSSW_5_3_32/MG5_aMC_v2_6_7/MY_FIRST_MG5_RUN/Cards/run_card.dat

Do you want to edit a card (press enter to bypass editing)?

/------------------------------------------------------------\

|  1. param : param_card.dat                                 |

|  2. run   : run_card.dat                                   |

\------------------------------------------------------------/

 you can also

   - enter the path to a valid card or banner.

   - use the 'set' command to modify a parameter directly.

     The set option works only for param_card and run_card.

     Type 'help set' for more information on this command.

   - call an external program (ASperGE/MadWidth/...).

     Type 'help' for the list of available command

 [0, done, 1, param, 2, run, enter path] 

>0

INFO: Update the dependent parameter of the param_card.dat 

Generating 10000 events with run name run_01

survey  run_01 

INFO: compile directory 

compile Source Directory

Using random number seed offset = 21

INFO: Running Survey 

Creating Jobs

Working on SubProcesses

INFO:     P1_gg_ttx  

INFO:     P1_qq_ttx  

INFO:     P2_gq_wpq_wp_lvl  

INFO:     P2_qq_wpg_wp_lvl  

INFO:  Idle: 1,  Running: 0,  Completed: 3 [ current time: 12h22 ] 

INFO:  Idle: 0,  Running: 1,  Completed: 3 [ current time: 12h22 ] 

INFO:  Idle: 0,  Running: 0,  Completed: 4 [  0.92s  ] 

INFO: End survey 

refine 10000

Creating Jobs

INFO: Refine results to 10000 

INFO: Generating 10000.0 unweigthed events. 

INFO: Effective Luminosity 2.33720491327 pb^-1 

INFO: need to improve 5 channels 

Current estimate of cross-section: 5134.338 +- 75.8652574554

    P1_gg_ttx 

    P1_qq_ttx 

    P2_gq_wpq_wp_lvl 

    P2_qq_wpg_wp_lvl 

INFO:  Idle: 13,  Running: 1,  Completed: 0 [ current time: 12h22 ] 

INFO:  Idle: 11,  Running: 1,  Completed: 2 [  2.5s  ] 

INFO:  Idle: 9,  Running: 1,  Completed: 4 [  14s  ] 

INFO:  Idle: 7,  Running: 1,  Completed: 6 [  17.8s  ] 

INFO:  Idle: 5,  Running: 1,  Completed: 8 [  21.6s  ] 

INFO:  Idle: 3,  Running: 1,  Completed: 10 [  25.3s  ] 

INFO:  Idle: 0,  Running: 0,  Completed: 14 [  28s  ] 

INFO: Combining runs 

INFO: finish refine 

refine 10000

Creating Jobs

INFO: Refine results to 10000 

INFO: Generating 10000.0 unweigthed events. 

INFO: Effective Luminosity 2.3200384817 pb^-1 

INFO: need to improve 0 channels 

Current estimate of cross-section: 5172.328 +- 15.5613446801

    P1_gg_ttx 

    P1_qq_ttx 

    P2_gq_wpq_wp_lvl 

    P2_qq_wpg_wp_lvl 

INFO:  Idle: 0,  Running: 0,  Completed: 0 [ current time: 12h23 ] 

INFO: Combining runs 

INFO: finish refine 

INFO: Combining Events 

  === Results Summary for run: run_01 tag: tag_1 ===

     Cross-section :   5172 +- 15.56 pb

     Nb of events :  10000

INFO: No version of lhapdf. Can not run systematics computation 

store_events

INFO: Storing parton level results 

INFO: End Parton 

reweight -from_cards

decay_events -from_cards

INFO: storing files of previous run 

INFO: Done 

quit

INFO:  

more information in /home/cms-opendata/CMSSW_5_3_32/MG5_aMC_v2_6_7/MY_FIRST_MG5_RUN/index.html

MG5_aMC> 

WebPage example

gunzip unweighted_events.lhe.gz

python lhe2root.py unweighted_events.lhe MY_FIRST_MG5_RUN.root

 root -l MY_FIRST_MG5_RUN.root 

TBrowser t #to open in root

Generated graphs within TBrowser t#include "Riostream.h" #include <iostream> #include <fstream> /* some code to help us understand why and how resolutions affect our ability to measure the mass of a particle */ // main routine void resolutions(){ // set debug mode bool idebug=false; // set up the display gROOT->Reset(); gROOT->SetStyle("Plain"); gStyle->SetOptStat(0); gStyle->SetOptFit(11); gStyle->SetErrorX(0); TCanvas *c1 = new TCanvas("c1","Our data",200,10,700,500); c1->SetFillColor(0); c1->SetBorderMode(0); c1->SetBorderSize(1); c1->SetFrameBorderMode(0); TCanvas *c2 = new TCanvas("c2","Our data",200,10,700,500); c2->SetFillColor(0); c2->SetBorderMode(0); c2->SetBorderSize(1); c2->SetFrameBorderMode(0); // book some histograms TH1F *histo1 = new TH1F("histo1","mass",1000,0.,200.); TH1F *histo2 = new TH1F("histo2","transverse momentum",1000,0.,200.); // setup random number generator gRandom->SetSeed(); // get the secret parameters ifstream myfile; myfile.open("secretparameters.txt"); double mass, resE, resA; myfile>>mass>>resE>>resA; myfile.close(); // make N bosons at rest // each boson will decay to 2 back-to-back daughter particles // they will be back-to-back by conservation of momentum, since the momentum of the mother was zero // assume that the mass of the daughter particles is small compared to the mother mass so we can // assume that their energy will be large compared to their mass and we can treat them as massless. // and thus their energy is equal to the magnitude of their momenta // by conservations of energy, their energy // work in 2D. Can always choose my coordinate system so that the 2 daughters are in the same plane int N=100; double etrue1,etrue2,phitrue1,phitrue2; // true energy of the 2 daughters double e1,px1,py1,phi1; // smeared 4-momenta of daughter 1 double e2,px2,py2,phi2; // smeared 4-momenta of daughter 2 double masssmeared, p_T; for(int i=0;i<N;i++) { // set true energy etrue1=mass/2.; etrue2=mass/2.; if(idebug) cout<<"etrue "<<etrue1<<" "<<etrue2<<endl; // choose phi for daughter 1 and daughter 2 phitrue1=2*TMath::Pi()*gRandom->Rndm(); phitrue2 = phitrue1+TMath::Pi(); if(idebug) cout<<"phitrue "<<phitrue1<<" "<<phitrue2<<endl; // smear true energy with resolution of detector e1=etrue1+resE*gRandom->Gaus(0.,1.); e2=etrue2+resE*gRandom->Gaus(0.,1.); if(idebug) cout<<"e "<<e1<<" "<<e2<<endl; //smear angles with resolution of the detector phi1=phitrue1+resA*gRandom->Gaus(0.,1.); phi2=phitrue2+resA*gRandom->Gaus(0.,1.); if(idebug) cout<<"phi "<<phi1<<" "<<phi2<<endl; //calculate 4 momenta after smearing px1=e1*cos(phi1); py1=e1*sin(phi1); px2=e2*cos(phi2); py2=e2*sin(phi2); if(idebug) cout<<"pxs "<<px1<<" "<<py1<<" "<<px2<<" "<<py2<<endl; // calculate smeared mass masssmeared=sqrt((e1+e2)*(e1+e2) - (px1+px2)*(px1+px2) - (py1+py2)*(py1+py2)); //my terrible code p_T = sqrt( (px1+px2)*(px1+px2)+(py1+py2)*(py1+py2) ); histo2->Fill(p_T); if(idebug) cout<<"masssmeared "<<masssmeared<<endl; histo1->Fill(masssmeared); } c1->cd(); histo1->Draw(""); c1->Update(); c2->cd(); histo2->Draw(""); c2->Update(); c1->SaveAs("c1.gif"); c2->SaveAs("c2.gif"); }

HW 8 

Tlorentz code

#define HiggsAnalysis_cxx

#include "HiggsAnalysis.h"

#include <TH2.h>

#include <TStyle.h>

#include <TCanvas.h>

void HiggsAnalysis::Loop()

{

// In a ROOT session, you can do:

// Root > .L HiggsAnalysis.C

// Root > HiggsAnalysis t

// Root > t.GetEntry(12); // Fill t data members with entry number 12

// Root > t.Show(); // Show values of entry 12

// Root > t.Show(16); // Read and show values of entry 16

// Root > t.Loop(); // Loop on all entries

//

// This is the loop skeleton where:

// jentry is the global entry number in the chain

// ientry is the entry number in the current Tree

// Note that the argument to GetEntry must be:

// jentry for TChain::GetEntry

// ientry for TTree::GetEntry and TBranch::GetEntry

//

// To read only selected branches, Insert statements like:

// METHOD1:

// fChain->SetBranchStatus("*",0); // disable all branches

// fChain->SetBranchStatus("branchname",1); // activate branchname

// METHOD2: replace line

// fChain->GetEntry(jentry); //read all branches

//by b_branchname->GetEntry(ientry); //read only this branch

if (fChain == 0) return;

Long64_t nentries = fChain->GetEntriesFast();

Long64_t nbytes = 0, nb = 0;

TFile* output = TFile::Open("Dielectron_MC.root", "RECREATE"); // "RECREATE" would produce a new root file with name Dielectron_MC.root every time you run the code

TH1F* Z_ee = new TH1F("Z_ee", "Di-electron candidate invariant mass", 200, 0, 200);

double el1mt = 0.0;

double el1pt = 0.0;

double el1eta = 0.0;

for (Long64_t jentry=0; jentry<nentries;jentry++) {

Long64_t ientry = LoadTree(jentry);

if (ientry < 0) break;

nb = fChain->GetEntry(jentry); nbytes += nb;

// if (Cut(ientry) < 0) continue;

TLorentzVector el1, el2;

el1.SetPtEtaPhiM(f_lept1_pt, f_lept1_eta, f_lept1_phi, 0.0);

el2.SetPtEtaPhiM(f_lept2_pt, f_lept2_eta, f_lept2_phi, 0.0);

TLorentzVector zCandidate = el1 + el2;

Z_ee->Fill(zCandidate.M());

el1mt = el1.Mt();

cout << el1mt << endl;

}

Z_ee->Write();

output->Close();

}

LOG 9

/**********************************************************************************

* Project : TMVA - a ROOT-integrated toolkit for multivariate data analysis *

* Package : TMVA *

* Root Macro: TMVAClassification *

* *

* This macro provides examples for the training and testing of the *

* TMVA classifiers. *

* *

* As input data is used a toy-MC sample consisting of four Gaussian-distributed *

11 * and linearly correlated input variables. *

12 * *

13 * The methods to be used can be switched on and off by means of booleans, or *

14 * via the prompt command, for example: *

15 * *

16 * root -l ./TMVAClassification.C\(\"my method\"\) 

// or to use default select method root -l ./TMVAClassification.C(\"\"\) 

17 * *

18 * (note that the backslashes are mandatory) *

*

21 * The output file "TMVA.root" can be analysed with the use of dedicated *

22 * macros (simply say: root -l <macro.C>), which can be conveniently *

23 * invoked through a GUI that will appear at the end of the run of this macro. *

24 * Launch the GUI via the command: *

25 * *

26 * root -l ./TMVAGui.C *

27 * * *

9 * If no method given, a default set is of classifiers is used *

0 **********************************************************************************/ 

#include <cstdlib>

#include <iostream>

#include <map>

#include <string>

#include "TChain.h"

#include "TFile.h"

#include "TTree.h"

#include "TString.h"

#include "TObjString.h"

#include "TSystem.h"

#include "TROOT.h"

#include "TMVA/Factory.h"

#include "TMVA/Tools.h"

//#include "TMVA/TMVAGui.h"

int TMVAClassification( TString myMethodList = "BDTB, MLP, Likelihood" )

{

// The explicit loading of the shared libTMVA is done in TMVAlogon.C, defined in .rootrc

// if you use your private .rootrc, or run from a different directory, please copy the

// corresponding lines from .rootrc

// methods to be processed can be given as an argument; use format:

//

// mylinux~> root -l TMVAClassification.CUndefined control sequence \"

//

// if you like to use a method via the plugin mechanism, we recommend using

//

// mylinux~> root -l TMVAClassification.CUndefined control sequence \"

// (an example is given for using the BDT as plugin (see below),

// but of course the real application is when you write your own

// method based)

//---------------------------------------------------------------

// This loads the library

TMVA::Tools::Instance();

// Default MVA methods to be trained + tested

std::map<std::string,int> Use;

// --- Cut optimisation

Use["Cuts"] = 1;

Use["CutsD"] = 1;

Use["CutsPCA"] = 0;

Use["CutsGA"] = 0;

Use["CutsSA"] = 0;

// 

// --- 1-dimensional likelihood ("naive Bayes estimator")

Use["Likelihood"] = 1;

Use["LikelihoodD"] = 0; // the "D" extension indicates decorrelated input variables (see option strings)

Use["LikelihoodPCA"] = 0; // the "PCA" extension indicates PCA-transformed input variables (see option strings)

Use["LikelihoodKDE"] = 0;

Use["LikelihoodMIX"] = 0;

//

// --- Mutidimensional likelihood and Nearest-Neighbour methods

Use["PDERS"] = 0;

Use["PDERSD"] = 0;

Use["PDERSPCA"] = 0;

Use["PDEFoam"] = 0;

Use["PDEFoamBoost"] = 0; // uses generalised MVA method boosting

Use["KNN"] = 1; // k-nearest neighbour method

//

// --- Linear Discriminant Analysis

Use["LD"] = 0; // Linear Discriminant identical to Fisher

Use["Fisher"] = 0;

Use["FisherG"] = 0;

Use["BoostedFisher"] = 0; // uses generalised MVA method boosting

Use["HMatrix"] = 0;

//

// --- Function Discriminant analysis

Use["FDA_GA"] = 0; // minimisation of user-defined function using Genetics Algorithm

Use["FDA_SA"] = 0;

Use["FDA_MC"] = 0;

Use["FDA_MT"] = 0;

Use["FDA_GAMT"] = 0;

Use["FDA_MCMT"] = 0;

//

// --- Neural Networks (all are feed-forward Multilayer Perceptrons)

Use["MLP"] = 1; // Recommended ANN

Use["MLPBFGS"] = 0; // Recommended ANN with optional training method

Use["MLPBNN"] = 0; // Recommended ANN with BFGS training method and bayesian regulator

Use["CFMlpANN"] = 0; // Depreciated ANN from ALEPH

Use["TMlpANN"] = 0; // ROOT's own ANN

//

// --- Support Vector Machine 

Use["SVM"] = 1;

// 

// --- Boosted Decision Trees

Use["BDT"] = 1; // uses Adaptive Boost

Use["BDTG"] = 0; // uses Gradient Boost

Use["BDTB"] = 0; // uses Bagging

Use["BDTD"] = 0; // decorrelation + Adaptive Boost

Use["BDTF"] = 0; // allow usage of fisher discriminant for node splitting 

// 

// --- Friedman's RuleFit method, ie, an optimised series of cuts ("rules")

Use["RuleFit"] = 1;

// ---------------------------------------------------------------

std::cout << std::endl;

std::cout << "==> Start TMVAClassification" << std::endl;

// Select methods (don't look at this code - not of interest)

if (myMethodList != "") {

for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0;

std::vector<TString> mlist = TMVA::gTools().SplitString( myMethodList, ',' );

for (UInt_t i=0; i<mlist.size(); i++) {

std::string regMethod(mlist[i]);

if (Use.find(regMethod) == Use.end()) {

std::cout << "Method \"" << regMethod << "\" not known in TMVA under this name. Choose among the following:" << std::endl;

for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) std::cout << it->first << " ";

std::cout << std::endl;

return 1;

}

Use[regMethod] = 1;

}

}

// --------------------------------------------------------------------------------------------------

// --- Here the preparation phase begins

// Create a ROOT output file where TMVA will store ntuples, histograms, etc.

TString outfileName( "TMVA.root" );

TFile* outputFile = TFile::Open( outfileName, "RECREATE" );

// Create the factory object. Later you can choose the methods

// whose performance you'd like to investigate. The factory is 

// the only TMVA object you have to interact with

//

// The first argument is the base of the name of all the

// weightfiles in the directory weight/

//

// The second argument is the output file for the training results

// All TMVA output can be suppressed by removing the "!" (not) in

// front of the "Silent" argument in the option string

TMVA::Factory *factory = new TMVA::Factory( "TMVAClassification", outputFile,

"!V:!Silent:Color:DrawProgressBar:Transformations=I;D;P;G,D:AnalysisType=Classification" );

// If you wish to modify default settings

// (please check "src/Config.h" to see all available global options)

// (TMVA::gConfig().GetVariablePlotting()).fTimesRMS = 8.0;

// (TMVA::gConfig().GetIONames()).fWeightFileDir = "myWeightDirectory";

// Define the input variables that shall be used for the MVA training

// note that you may also use variable expressions, such as: "3*var1/var2*abs(var3)"

// [all types of expressions that can also be parsed by TTree::Draw( "expression" )]

factory->AddVariable( "myvar1 := var1+var2", 'F' );

factory->AddVariable( "myvar2 := var1-var2", "Expression 2", "", 'F' );

factory->AddVariable( "var3", "Variable 3", "units", 'F' );

factory->AddVariable( "var4", "Variable 4", "units", 'F' );

// You can add so-called "Spectator variables", which are not used in the MVA training,

// but will appear in the final "TestTree" produced by TMVA. This TestTree will contain the

// input variables, the response values of all trained MVAs, and the spectator variables

factory->AddSpectator( "spec1 := var1*2", "Spectator 1", "units", 'F' );

factory->AddSpectator( "spec2 := var1*3", "Spectator 2", "units", 'F' );

// Read training and test data

// (it is also possible to use ASCII format as input -> see TMVA Users Guide)

TString fname = "./tmva_class_example.root";

if (gSystem->AccessPathName( fname )) // file does not exist in local directory

gSystem->Exec("curl -O http://root.cern.ch/files/tmva_class_example.root");

TFile *input = TFile::Open( fname );

std::cout << "--- TMVAClassification : Using input file: " << input->GetName() << std::endl;

// --- Register the training and test trees

TTree *signal = (TTree*)input->Get("TreeS");

TTree *background = (TTree*)input->Get("TreeB");

// global event weights per tree (see below for setting event-wise weights)

Double_t signalWeight = 1.0;

Double_t backgroundWeight = 1.0;

// You can add an arbitrary number of signal or background trees

factory->AddSignalTree ( signal, signalWeight );

factory->AddBackgroundTree( background, backgroundWeight );

// To give different trees for training and testing, do as follows:

// factory->AddSignalTree( signalTrainingTree, signalTrainWeight, "Training" );

// factory->AddSignalTree( signalTestTree, signalTestWeight, "Test" );

// Use the following code instead of the above two or four lines to add signal and background

// training and test events "by hand"

// NOTE that in this case one should not give expressions (such as "var1+var2") in the input

// variable definition, but simply compute the expression before adding the event

//

// // --- begin ----------------------------------------------------------

// std::vector<Double_t> vars( 4 ); // vector has size of number of input variables

// Float_t treevars[4], weight;

// 

// // Signal

// for (UInt_t ivar=0; ivar<4; ivar++) signal->SetBranchAddress( Form( "var%i", ivar+1 ), &(treevars[ivar]) );

// for (UInt_t i=0; i<signal->GetEntries(); i++) {

// signal->GetEntry(i);

// for (UInt_t ivar=0; ivar<4; ivar++) vars[ivar] = treevars[ivar];

// // add training and test events; here: first half is training, second is testing

// // note that the weight can also be event-wise

// if (i < signal->GetEntries()/2.0) factory->AddSignalTrainingEvent( vars, signalWeight );

// else factory->AddSignalTestEvent ( vars, signalWeight );

// }

// 

// // Background (has event weights)

// background->SetBranchAddress( "weight", &weight );

// for (UInt_t ivar=0; ivar<4; ivar++) background->SetBranchAddress( Form( "var%i", ivar+1 ), &(treevars[ivar]) );

// for (UInt_t i=0; i<background->GetEntries(); i++) {

// background->GetEntry(i);

// for (UInt_t ivar=0; ivar<4; ivar++) vars[ivar] = treevars[ivar];

// // add training and test events; here: first half is training, second is testing

// // note that the weight can also be event-wise

// if (i < background->GetEntries()/2) factory->AddBackgroundTrainingEvent( vars, backgroundWeight*weight );

// else factory->AddBackgroundTestEvent ( vars, backgroundWeight*weight );

// }

// --- end ------------------------------------------------------------

//

// --- end of tree registration 

// Set individual event weights (the variables must exist in the original TTree)

// for signal : factory->SetSignalWeightExpression ("weight1*weight2");

// for background: factory->SetBackgroundWeightExpression("weight1*weight2");

factory->SetBackgroundWeightExpression( "weight" );

// Apply additional cuts on the signal and background samples (can be different)

TCut mycuts = ""; // for example: TCut mycuts = "abs(var1)<0.5 && abs(var2-0.5)<1";

TCut mycutb = ""; // for example: TCut mycutb = "abs(var1)<0.5";

// Tell the factory how to use the training and testing events

//

// If no numbers of events are given, half of the events in the tree are used 

// for training, and the other half for testing:

// factory->PrepareTrainingAndTestTree( mycut, "SplitMode=random:!V" );

// To also specify the number of testing events, use:

// factory->PrepareTrainingAndTestTree( mycut,

// "NSigTrain=3000:NBkgTrain=3000:NSigTest=3000:NBkgTest=3000:SplitMode=Random:!V" );

factory->PrepareTrainingAndTestTree( mycuts, mycutb,

"nTrain_Signal=1000:nTrain_Background=1000:SplitMode=Random:NormMode=NumEvents:!V" );

// ---- Book MVA methods

//

// Please lookup the various method configuration options in the corresponding cxx files, eg:

// src/MethoCuts.cxx, etc, or here: http://tmva.sourceforge.net/optionRef.html

// it is possible to preset ranges in the option string in which the cut optimisation should be done:

// "...:CutRangeMin[2]=-1:CutRangeMax[2]=1"...", where [2] is the third input variable

// Cut optimisation

if (Use["Cuts"])

factory->BookMethod( TMVA::Types::kCuts, "Cuts",

"!H:!V:FitMethod=MC:EffSel:SampleSize=200000:VarProp=FSmart" );

if (Use["CutsD"])

factory->BookMethod( TMVA::Types::kCuts, "CutsD",

"!H:!V:FitMethod=MC:EffSel:SampleSize=200000:VarProp=FSmart:VarTransform=Decorrelate" );

if (Use["CutsPCA"])

factory->BookMethod( TMVA::Types::kCuts, "CutsPCA",

"!H:!V:FitMethod=MC:EffSel:SampleSize=200000:VarProp=FSmart:VarTransform=PCA" );

if (Use["CutsGA"])

factory->BookMethod( TMVA::Types::kCuts, "CutsGA",

"H:!V:FitMethod=GA:CutRangeMin[0]=-10:CutRangeMax[0]=10:VarProp[1]=FMax:EffSel:Steps=30:Cycles=3:PopSize=400:SC_steps=10:SC_rate=5:SC_factor=0.95" );

if (Use["CutsSA"])

factory->BookMethod( TMVA::Types::kCuts, "CutsSA",

"!H:!V:FitMethod=SA:EffSel:MaxCalls=150000:KernelTemp=IncAdaptive:InitialTemp=1e+6:MinTemp=1e-6:Eps=1e-10:UseDefaultScale" );

// Likelihood ("naive Bayes estimator")

if (Use["Likelihood"])

factory->BookMethod( TMVA::Types::kLikelihood, "Likelihood",

"H:!V:TransformOutput:PDFInterpol=Spline2:NSmoothSig[0]=20:NSmoothBkg[0]=20:NSmoothBkg[1]=10:NSmooth=1:NAvEvtPerBin=50" );

// Decorrelated likelihood

if (Use["LikelihoodD"])

factory->BookMethod( TMVA::Types::kLikelihood, "LikelihoodD",

"!H:!V:TransformOutput:PDFInterpol=Spline2:NSmoothSig[0]=20:NSmoothBkg[0]=20:NSmooth=5:NAvEvtPerBin=50:VarTransform=Decorrelate" );

// PCA-transformed likelihood

if (Use["LikelihoodPCA"])

factory->BookMethod( TMVA::Types::kLikelihood, "LikelihoodPCA",

"!H:!V:!TransformOutput:PDFInterpol=Spline2:NSmoothSig[0]=20:NSmoothBkg[0]=20:NSmooth=5:NAvEvtPerBin=50:VarTransform=PCA" ); 

// Use a kernel density estimator to approximate the PDFs

if (Use["LikelihoodKDE"])

factory->BookMethod( TMVA::Types::kLikelihood, "LikelihoodKDE",

"!H:!V:!TransformOutput:PDFInterpol=KDE:KDEtype=Gauss:KDEiter=Adaptive:KDEFineFactor=0.3:KDEborder=None:NAvEvtPerBin=50" ); 

// Use a variable-dependent mix of splines and kernel density estimator

if (Use["LikelihoodMIX"])

factory->BookMethod( TMVA::Types::kLikelihood, "LikelihoodMIX",

"!H:!V:!TransformOutput:PDFInterpolSig[0]=KDE:PDFInterpolBkg[0]=KDE:PDFInterpolSig[1]=KDE:PDFInterpolBkg[1]=KDE:PDFInterpolSig[2]=Spline2:PDFInterpolBkg[2]=Spline2:PDFInterpolSig[3]=Spline2:PDFInterpolBkg[3]=Spline2:KDEtype=Gauss:KDEiter=Nonadaptive:KDEborder=None:NAvEvtPerBin=50" ); 

// Test the multi-dimensional probability density estimator

// here are the options strings for the MinMax and RMS methods, respectively:

// "!H:!V:VolumeRangeMode=MinMax:DeltaFrac=0.2:KernelEstimator=Gauss:GaussSigma=0.3" );

// "!H:!V:VolumeRangeMode=RMS:DeltaFrac=3:KernelEstimator=Gauss:GaussSigma=0.3" );

if (Use["PDERS"])

factory->BookMethod( TMVA::Types::kPDERS, "PDERS",

"!H:!V:NormTree=T:VolumeRangeMode=Adaptive:KernelEstimator=Gauss:GaussSigma=0.3:NEventsMin=400:NEventsMax=600" );

if (Use["PDERSD"])

factory->BookMethod( TMVA::Types::kPDERS, "PDERSD",

"!H:!V:VolumeRangeMode=Adaptive:KernelEstimator=Gauss:GaussSigma=0.3:NEventsMin=400:NEventsMax=600:VarTransform=Decorrelate" );

if (Use["PDERSPCA"])

factory->BookMethod( TMVA::Types::kPDERS, "PDERSPCA",

"!H:!V:VolumeRangeMode=Adaptive:KernelEstimator=Gauss:GaussSigma=0.3:NEventsMin=400:NEventsMax=600:VarTransform=PCA" );

// Multi-dimensional likelihood estimator using self-adapting phase-space binning

if (Use["PDEFoam"])

factory->BookMethod( TMVA::Types::kPDEFoam, "PDEFoam",

"!H:!V:SigBgSeparate=F:TailCut=0.001:VolFrac=0.0666:nActiveCells=500:nSampl=2000:nBin=5:Nmin=100:Kernel=None:Compress=T" );

if (Use["PDEFoamBoost"])

factory->BookMethod( TMVA::Types::kPDEFoam, "PDEFoamBoost",

"!H:!V:Boost_Num=30:Boost_Transform=linear:SigBgSeparate=F:MaxDepth=4:UseYesNoCell=T:DTLogic=MisClassificationError:FillFoamWithOrigWeights=F:TailCut=0:nActiveCells=500:nBin=20:Nmin=400:Kernel=None:Compress=T" );

// K-Nearest Neighbour classifier (KNN)

if (Use["KNN"])

factory->BookMethod( TMVA::Types::kKNN, "KNN",

"H:nkNN=20:ScaleFrac=0.8:SigmaFact=1.0:Kernel=Gaus:UseKernel=F:UseWeight=T:!Trim" );

// H-Matrix (chi2-squared) method

if (Use["HMatrix"])

factory->BookMethod( TMVA::Types::kHMatrix, "HMatrix", "!H:!V:VarTransform=None" );

// Linear discriminant (same as Fisher discriminant)

if (Use["LD"])

factory->BookMethod( TMVA::Types::kLD, "LD", "H:!V:VarTransform=None:CreateMVAPdfs:PDFInterpolMVAPdf=Spline2:NbinsMVAPdf=50:NsmoothMVAPdf=10" );

// Fisher discriminant (same as LD)

if (Use["Fisher"])

factory->BookMethod( TMVA::Types::kFisher, "Fisher", "H:!V:Fisher:VarTransform=None:CreateMVAPdfs:PDFInterpolMVAPdf=Spline2:NbinsMVAPdf=50:NsmoothMVAPdf=10" );

// Fisher with Gauss-transformed input variables

if (Use["FisherG"])

factory->BookMethod( TMVA::Types::kFisher, "FisherG", "H:!V:VarTransform=Gauss" );

// Composite classifier: ensemble (tree) of boosted Fisher classifiers

if (Use["BoostedFisher"])

factory->BookMethod( TMVA::Types::kFisher, "BoostedFisher", 

"H:!V:Boost_Num=20:Boost_Transform=log:Boost_Type=AdaBoost:Boost_AdaBoostBeta=0.2:!Boost_DetailedMonitoring" );

// Function discrimination analysis (FDA) -- test of various fitters - the recommended one is Minuit (or GA or SA)

if (Use["FDA_MC"])

factory->BookMethod( TMVA::Types::kFDA, "FDA_MC",

"H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=MC:SampleSize=100000:Sigma=0.1" );

if (Use["FDA_GA"]) // can also use Simulated Annealing (SA) algorithm (see Cuts_SA options])

factory->BookMethod( TMVA::Types::kFDA, "FDA_GA",

"H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=GA:PopSize=300:Cycles=3:Steps=20:Trim=True:SaveBestGen=1" );

if (Use["FDA_SA"]) // can also use Simulated Annealing (SA) algorithm (see Cuts_SA options])

factory->BookMethod( TMVA::Types::kFDA, "FDA_SA",

"H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=SA:MaxCalls=15000:KernelTemp=IncAdaptive:InitialTemp=1e+6:MinTemp=1e-6:Eps=1e-10:UseDefaultScale" );

if (Use["FDA_MT"])

factory->BookMethod( TMVA::Types::kFDA, "FDA_MT",

"H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=MINUIT:ErrorLevel=1:PrintLevel=-1:FitStrategy=2:UseImprove:UseMinos:SetBatch" );

if (Use["FDA_GAMT"])

factory->BookMethod( TMVA::Types::kFDA, "FDA_GAMT",

"H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=GA:Converger=MINUIT:ErrorLevel=1:PrintLevel=-1:FitStrategy=0:!UseImprove:!UseMinos:SetBatch:Cycles=1:PopSize=5:Steps=5:Trim" );

if (Use["FDA_MCMT"])

factory->BookMethod( TMVA::Types::kFDA, "FDA_MCMT",

"H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=MC:Converger=MINUIT:ErrorLevel=1:PrintLevel=-1:FitStrategy=0:!UseImprove:!UseMinos:SetBatch:SampleSize=20" );

// TMVA ANN: MLP (recommended ANN) -- all ANNs in TMVA are Multilayer Perceptrons

if (Use["MLP"])

factory->BookMethod( TMVA::Types::kMLP, "MLP", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:!UseRegulator" );

if (Use["MLPBFGS"])

factory->BookMethod( TMVA::Types::kMLP, "MLPBFGS", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:TrainingMethod=BFGS:!UseRegulator" );

if (Use["MLPBNN"])

factory->BookMethod( TMVA::Types::kMLP, "MLPBNN", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:TrainingMethod=BFGS:UseRegulator" ); // BFGS training with bayesian regulators

// CF(Clermont-Ferrand)ANN

if (Use["CFMlpANN"])

factory->BookMethod( TMVA::Types::kCFMlpANN, "CFMlpANN", "!H:!V:NCycles=2000:HiddenLayers=N+1,N" ); // n_cycles:#nodes:#nodes:... 

// Tmlp(Root)ANN

if (Use["TMlpANN"])

factory->BookMethod( TMVA::Types::kTMlpANN, "TMlpANN", "!H:!V:NCycles=200:HiddenLayers=N+1,N:LearningMethod=BFGS:ValidationFraction=0.3" ); // n_cycles:#nodes:#nodes:...

// Support Vector Machine

if (Use["SVM"])

factory->BookMethod( TMVA::Types::kSVM, "SVM", "Gamma=0.25:Tol=0.001:VarTransform=Norm" );

// Boosted Decision Trees

if (Use["BDTG"]) // Gradient Boost

factory->BookMethod( TMVA::Types::kBDT, "BDTG",

"!H:!V:NTrees=1000:MinNodeSize=2.5%:BoostType=Grad:Shrinkage=0.10:UseBaggedBoost:BaggedSampleFraction=0.5:nCuts=20:MaxDepth=2" );

if (Use["BDT"]) // Adaptive Boost

//factory->BookMethod( TMVA::Types::kBDT, "BDT",

// "!H:!V:NTrees=850:MinNodeSize=2.5%:MaxDepth=3:BoostType=AdaBoost:AdaBoostBeta=0.5:UseBaggedBoost:BaggedSampleFraction=0.5:SeparationType=GiniIndex:nCuts=20" );

factory->BookMethod( TMVA::Types::kBDT, "BDT",

"!H:!V:NTrees=850:MaxDepth=3:SeparationType=GiniIndex:nCuts=20" );

if (Use["BDTB"]) // Bagging

factory->BookMethod( TMVA::Types::kBDT, "BDTB",

"!H:!V:NTrees=400:BoostType=Bagging:SeparationType=GiniIndex:nCuts=20" );

if (Use["BDTD"]) // Decorrelation + Adaptive Boost

factory->BookMethod( TMVA::Types::kBDT, "BDTD",

"!H:!V:NTrees=400:MinNodeSize=5%:MaxDepth=3:BoostType=AdaBoost:SeparationType=GiniIndex:nCuts=20:VarTransform=Decorrelate" );

if (Use["BDTF"]) // Allow Using Fisher discriminant in node splitting for (strong) linearly correlated variables

factory->BookMethod( TMVA::Types::kBDT, "BDTMitFisher",

"!H:!V:NTrees=50:MinNodeSize=2.5%:UseFisherCuts:MaxDepth=3:BoostType=AdaBoost:AdaBoostBeta=0.5:SeparationType=GiniIndex:nCuts=20" );

// RuleFit -- TMVA implementation of Friedman's method

if (Use["RuleFit"])

factory->BookMethod( TMVA::Types::kRuleFit, "RuleFit",

"H:!V:RuleFitModule=RFTMVA:Model=ModRuleLinear:MinImp=0.001:RuleMinDist=0.001:NTrees=20:fEventsMin=0.01:fEventsMax=0.5:GDTau=-1.0:GDTauPrec=0.01:GDStep=0.01:GDNSteps=10000:GDErrScale=1.02" );

// For an example of the category classifier usage, see: TMVAClassificationCategory

// --------------------------------------------------------------------------------------------------

// ---- Now you can optimize the setting (configuration) of the MVAs using the set of training events

// ---- STILL EXPERIMENTAL and only implemented for BDT's ! 

// factory->OptimizeAllMethods("SigEffAt001","Scan");

// factory->OptimizeAllMethods("ROCIntegral","FitGA");

// --------------------------------------------------------------------------------------------------

// ---- Now you can tell the factory to train, test, and evaluate the MVAs

// Train MVAs using the set of training events

factory->TrainAllMethods();

// ---- Evaluate all MVAs using the set of test events

factory->TestAllMethods();

// ----- Evaluate and compare performance of all configured MVAs

factory->EvaluateAllMethods();

// --------------------------------------------------------------

// Save the output

outputFile->Close();

std::cout << "==> Wrote root file: " << outputFile->GetName() << std::endl;

std::cout << "==> TMVAClassification is done!" << std::endl;

delete factory;

// Launch the GUI for the root macros

//if (!gROOT->IsBatch()) TMVA::TMVAGui( outfileName );

return 0;

}

int main( int argc, char** argv )

{

// Select methods (don't look at this code - not of interest)

TString methodList; 

for (int i=1; i<argc; i++) {

TString regMethod(argv[i]);

if(regMethod=="-b" || regMethod=="--batch") continue;

if (!methodList.IsNull()) methodList += TString(","); 

methodList += regMethod;

}

return TMVAClassification(methodList); 

}