source: branches/Mars_McMismatchStudy/MMCComp/MMCComp.cpp@ 18100

#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 + "/gan*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;

}