Index: /branches/Mars_McMismatchStudy/MMCComp/MMCComp.cpp =================================================================== --- /branches/Mars_McMismatchStudy/MMCComp/MMCComp.cpp (revision 18090) +++ /branches/Mars_McMismatchStudy/MMCComp/MMCComp.cpp (revision 18090) @@ -0,0 +1,263 @@ +#include "MMCComp.h" + + +MMCComp::MMCComp() +{ +} + +void MMCComp::Print() const +{ + +/* + for( int i = 0 ; i < cDATA->GetEntries() ; i++ ) + { + cDATA->GetEntry(i); + cout << i << " "; + cout << time->GetTime() << " "; + cout << time->GetMjd() << " "; + cout << endl; + } +*/ + +} + + +void MMCComp::loadData( string sDataDir ) +{ + + string sDataFiles = sDataDir + "2*summary.root"; + cDATA = new TChain("Events"); + + cDATA->Add( sDataFiles.c_str() ); + cout << "Adding " << sDataFiles.c_str() << endl; + cout << "Number of Data Events = " << cDATA->GetEntries() << endl; + + timeDATA = new MTime(); + hillasDATA = new MHillas(); + imageparDATA = new MImagePar(); + newimageparDATA = new MNewImagePar();; + + cDATA->SetBranchStatus("MTime.*",1); + cDATA->SetBranchAddress("MTime.",&timeDATA); + + cDATA->SetBranchStatus("MHillas.*",1); + cDATA->SetBranchAddress("MHillas.",&hillasDATA); + + cDATA->SetBranchStatus("MImagePar.*",1); + cDATA->SetBranchAddress("MImagePar.",&imageparDATA); + + cDATA->SetBranchStatus("MNewImagePar.*",1); + cDATA->SetBranchAddress("MNewImagePar.",&newimageparDATA); + +} + +void MMCComp::loadMC( string sMCDir ) +{ + + string sMCFiles = sMCDir + "/ganymedMC_Full_20140716-analysis.root"; + cMC = new TChain("Events"); + + cMC->Add( sMCFiles.c_str() ); + cout << "Adding " << sMCFiles.c_str() << endl; + cout << "Number of MC Events = " << cMC->GetEntries() << endl; + + //mcEvent = new MMcEvtBasic(); + hillasMC = new MHillas(); + imageparMC = new MImagePar(); + newimageparMC = new MNewImagePar();; + + //cMC->SetBranchStatus("MMcEvtBasic.*",1); + //cMC->SetBranchAddress("MMcEvtBasic.",&mcEvent); + + cMC->SetBranchStatus("MHillas.*",1); + cMC->SetBranchAddress("MHillas.",&hillasMC); + + cMC->SetBranchStatus("MImagePar.*",1); + cMC->SetBranchAddress("MImagePar.",&imageparMC); + + cMC->SetBranchStatus("MNewImagePar.*",1); + cMC->SetBranchAddress("MNewImagePar.",&newimageparMC); + + +} + +void MMCComp::compareSize( int N_Size, double MIN_Size, double MAX_Size ) +{ + + double min,max; + + min = 9999999.0; + max = -9999999.0; + + TH1F *hSizeDATA = new TH1F("hSizeDATA","",N_Size,MIN_Size,MAX_Size); + TH1F *hSizeMC = new TH1F("hSizeMC","",N_Size,MIN_Size,MAX_Size); + TH1F *hSizeRATIO = new TH1F("hSizeRATIO","",N_Size,MIN_Size,MAX_Size); + + for( int i = 0 ; i < cDATA->GetEntries() ; i++ ) + { + cDATA->GetEntry(i); + hSizeDATA->Fill( log10(hillasDATA->GetSize()) ); + if( log10(hillasDATA->GetSize()) < min ) min = log10(hillasDATA->GetSize()); + if( log10(hillasDATA->GetSize()) > max ) max = log10(hillasDATA->GetSize()); + } + + for( int i = 0 ; i < cMC->GetEntries() ; i++ ) + { + cMC->GetEntry(i); + hSizeMC->Fill( log10(hillasMC->GetSize()) ); + } + + hSizeMC->Scale( hSizeDATA->GetEntries()/hSizeMC->GetEntries() ); + + for( int i = 0 ; i < hSizeDATA->GetNbinsX() ; i++ ) + if( hSizeMC->GetBinContent(i) > 0. && hSizeDATA->GetBinContent(i) > 0. ) + { + hSizeRATIO->Fill((double)hSizeDATA->GetBinCenter(i),hSizeDATA->GetBinContent(i)/hSizeMC->GetBinContent(i)); + hSizeRATIO->SetBinError((double)hSizeDATA->GetBinCenter(i),1./( 1./hSizeDATA->GetBinContent(i) + 1./hSizeMC->GetBinContent(i)) ); + } + + hSizeRATIO->SetMaximum(2.); + hSizeRATIO->SetMinimum(0.); + hSizeRATIO->SetMarkerStyle(7); + + hSizeDATA->SetMarkerStyle(7); + hSizeDATA->SetMarkerColor(1); + + hSizeMC->SetLineColor(2); + + TCanvas *canSize = new TCanvas("canSize","canSize",0,0,500,500); + canSize->Divide(1,2); + canSize->cd(1); + canSize->cd(1)->SetLogy(); + + hSizeDATA->Draw("e"); + hSizeMC->Draw("same"); + + canSize->cd(2); + hSizeRATIO->Draw("e"); + hSizeRATIO->Fit("pol1","","same"); + + cout << "Data Size min = " << min << " max = " << max << endl; + +} + +void MMCComp::compareLength( int N_Length, double MIN_Length, double MAX_Length ) +{ + + double min,max; + + min = 9999999.0; + max = -9999999.0; + + TH1F *hLengthDATA = new TH1F("hLengthDATA","",N_Length,MIN_Length,MAX_Length); + TH1F *hLengthMC = new TH1F("hLengthMC","",N_Length,MIN_Length,MAX_Length); + TH1F *hLengthRATIO = new TH1F("hLengthRATIO","",N_Length,MIN_Length,MAX_Length); + + for( int i = 0 ; i < cDATA->GetEntries() ; i++ ) + { + cDATA->GetEntry(i); + hLengthDATA->Fill( log10(hillasDATA->GetLength()) ); + if( log10(hillasDATA->GetLength()) < min ) min = log10(hillasDATA->GetLength()); + if( log10(hillasDATA->GetLength()) > max ) max = log10(hillasDATA->GetLength()); + } + + for( int i = 0 ; i < cMC->GetEntries() ; i++ ) + { + cMC->GetEntry(i); + hLengthMC->Fill( log10(hillasMC->GetLength()) ); + } + + hLengthMC->Scale( hLengthDATA->GetEntries()/hLengthMC->GetEntries() ); + + for( int i = 0 ; i < hLengthDATA->GetNbinsX() ; i++ ) + if( hLengthMC->GetBinContent(i) > 0. && hLengthDATA->GetBinContent(i) > 0. ) + { + hLengthRATIO->Fill((double)hLengthDATA->GetBinCenter(i),hLengthDATA->GetBinContent(i)/hLengthMC->GetBinContent(i)); + hLengthRATIO->SetBinError((double)hLengthDATA->GetBinCenter(i),1./( 1./hLengthDATA->GetBinContent(i) + 1./hLengthMC->GetBinContent(i)) ); + } + + hLengthRATIO->SetMaximum(2.); + hLengthRATIO->SetMinimum(0.); + hLengthRATIO->SetMarkerStyle(7); + + hLengthDATA->SetMarkerStyle(7); + hLengthDATA->SetMarkerColor(1); + + hLengthMC->SetLineColor(2); + + TCanvas *canLength = new TCanvas("canLength","canLength",0,0,500,500); + canLength->Divide(1,2); + canLength->cd(1); + canLength->cd(1)->SetLogy(); + + hLengthDATA->Draw("e"); + hLengthMC->Draw("same"); + + canLength->cd(2); + hLengthRATIO->Draw("e"); + hLengthRATIO->Fit("pol1","","same"); + + cout << "Data Length min = " << min << " max = " << max << endl; + +} + +void MMCComp::compareWidth( int N_Width, double MIN_Width, double MAX_Width ) +{ + + double min,max; + + min = 9999999.0; + max = -9999999.0; + + TH1F *hWidthDATA = new TH1F("hWidthDATA","",N_Width,MIN_Width,MAX_Width); + TH1F *hWidthMC = new TH1F("hWidthMC","",N_Width,MIN_Width,MAX_Width); + TH1F *hWidthRATIO = new TH1F("hWidthRATIO","",N_Width,MIN_Width,MAX_Width); + + for( int i = 0 ; i < cDATA->GetEntries() ; i++ ) + { + cDATA->GetEntry(i); + hWidthDATA->Fill( log10(hillasDATA->GetWidth()) ); + if( log10(hillasDATA->GetWidth()) < min ) min = log10(hillasDATA->GetWidth()); + if( log10(hillasDATA->GetWidth()) > max ) max = log10(hillasDATA->GetWidth()); + } + + for( int i = 0 ; i < cMC->GetEntries() ; i++ ) + { + cMC->GetEntry(i); + hWidthMC->Fill( log10(hillasMC->GetWidth()) ); + } + + hWidthMC->Scale( hWidthDATA->GetEntries()/hWidthMC->GetEntries() ); + + for( int i = 0 ; i < hWidthDATA->GetNbinsX() ; i++ ) + if( hWidthMC->GetBinContent(i) > 0. && hWidthDATA->GetBinContent(i) > 0. ) + { + hWidthRATIO->Fill((double)hWidthDATA->GetBinCenter(i),hWidthDATA->GetBinContent(i)/hWidthMC->GetBinContent(i)); + hWidthRATIO->SetBinError((double)hWidthDATA->GetBinCenter(i),1./( 1./hWidthDATA->GetBinContent(i) + 1./hWidthMC->GetBinContent(i)) ); + } + + hWidthRATIO->SetMaximum(2.); + hWidthRATIO->SetMinimum(0.); + hWidthRATIO->SetMarkerStyle(7); + + hWidthDATA->SetMarkerStyle(7); + hWidthDATA->SetMarkerColor(1); + + hWidthMC->SetLineColor(2); + + TCanvas *canWidth = new TCanvas("canWidth","canWidth",0,0,500,500); + canWidth->Divide(1,2); + canWidth->cd(1); + canWidth->cd(1)->SetLogy(); + + hWidthDATA->Draw("e"); + hWidthMC->Draw("same"); + + canWidth->cd(2); + hWidthRATIO->Draw("e"); + hWidthRATIO->Fit("pol1","","same"); + + cout << "Data Width min = " << min << " max = " << max << endl; + +} + Index: /branches/Mars_McMismatchStudy/MMCComp/MMCComp.h =================================================================== --- /branches/Mars_McMismatchStudy/MMCComp/MMCComp.h (revision 18090) +++ /branches/Mars_McMismatchStudy/MMCComp/MMCComp.h (revision 18090) @@ -0,0 +1,64 @@ +#ifndef __MMCOMP_H__ +#define __MMCOMP_H__ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "TCanvas.h" + +#include "../../mbase/MTime.h" +#include "../../mimage/MHillas.h" +#include "../../mimage/MImagePar.h" +#include "../../mimage/MNewImagePar.h" +//#include "../../mmc/MMcEvtBasic.h" + +//#include "/home/guest/ghughes/FACT/Mars_McMismatchStudy/libmars.so" + +//#define N_SIZE 100 +//#define MIN_SIZE 1. +//#define MAX_SIZE 5. + +using namespace std; + +class MMCComp +{ + +private: + TChain *cDATA; + TChain *cMC; + + MTime *timeDATA; + MHillas *hillasDATA; + MImagePar *imageparDATA; + MNewImagePar *newimageparDATA; + + MHillas *hillasMC; + MImagePar *imageparMC; + MNewImagePar *newimageparMC; + //MMcEvtBasic *mcEvent; + +public: + MMCComp(); + void Print() const; + void loadData( string sDataDir ); + void loadMC( string sDataDir ); + + void compareSize( int N_SIZE, double MIN_SIZE, double MAX_SIZE ); + void compareLength( int N_Length, double MIN_Length, double MAX_Length ); + void compareWidth( int N_Width, double MIN_Width, double MAX_Width ); + +}; + +#endif + Index: /branches/Mars_McMismatchStudy/MMCComp/Makefile =================================================================== --- /branches/Mars_McMismatchStudy/MMCComp/Makefile (revision 18090) +++ /branches/Mars_McMismatchStudy/MMCComp/Makefile (revision 18090) @@ -0,0 +1,34 @@ +##CXXFLAGS = -g -fPIC -I$(ROOTSYS)/include -L../../libmars.so +CXXFLAGS = -g -fPIC -I$(ROOTSYS)/include +LDFLAGS = -g -shared +LD = g++ -fPIC + +TTT = ../../libmars.so + +HDRS = MMCComp.h mclib.h + +SRCS = MMCComp.cpp mclib.cpp + +OBJS = MMCComp.o mclib.o + +ROOTLIB = $(shell root-config --glibs) + +PROGRAM = mclib.so + +all: $(PROGRAM) + +$(PROGRAM): $(OBJS) + @echo "Linking $(PROGRAM) ..." + @/bin/rm -f $(PROGRAM) + @$(LD) $(LDFLAGS) $(OBJS) $(ROOTLIB) $(TTT) -o $(PROGRAM) + @chmod 555 $(PROGRAM) + @echo "done" + +clean:; @rm -f $(OBJS) core + +### +MMCComp.o: MMCComp.h + +mclib.cpp: MMCComp.h + @echo "Generating dictionary ..." + @rootcint mclib.cpp -c MMCComp.h Index: /branches/Mars_McMismatchStudy/MMCComp/README =================================================================== --- /branches/Mars_McMismatchStudy/MMCComp/README (revision 18090) +++ /branches/Mars_McMismatchStudy/MMCComp/README (revision 18090) @@ -0,0 +1,18 @@ +sharedlibrary for comparing Data and Monte Carlo + +> make + +to generate the mclib.so + +> root +[] .L mclib.so + +[] MMCComp *a = new MMCComp() +[] a->loadData("$MARSSYS/OUTPUT/ganymed/") +[] a->loadMC("") +[] a->compareSize(100,1.2,5.) + + + + +