EvtNT3.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//
// Module: EvtGen/EvtNT3.hh
//
// Description: Class to handle NTuple for three body decays
// usage:
// Users need to provide the MC and data file in NTuple style. The decays are described as
//                        A->  x1    +   x2     +  x3
// angular distribution:  costheta1, costheta2 ,costheta3
// mass distribution:     m12, m13, and m23, here, m12 ==> mass of x1 and x2, etc.
// these variables are defined as double type, the tree name is mc and data, repectively
// Modification history:
//
//    PING RG     September 11, 2010         Module created
//
//------------------------------------------------------------------------
#include "EvtNT3.hh"
#include <fstream>
#include <iostream>
#include <stdlib.h>
#include <string>
 
void EvtNT3::init() {
  max     = 0;
  Ncos    = 20;
  chainMC = new TChain( "mc" );
  chainDT = new TChain( "data" );
 
  chainMC->SetDirectory( 0 );
  chainDT->SetDirectory( 0 );
 
  chainMC->Add( mcfile );
  chainDT->Add( datafile );
  //--MC
  chainMC->SetBranchAddress( "costheta1", &costheta1 );
  chainMC->SetBranchAddress( "costheta2", &costheta2 );
  chainMC->SetBranchAddress( "costheta3", &costheta3 );
  chainMC->SetBranchAddress( "m12", &m12 );
  chainMC->SetBranchAddress( "m13", &m13 );
  chainMC->SetBranchAddress( "m23", &m23 );
  //--Data
  chainDT->SetBranchAddress( "costheta1", &costheta1 );
  chainDT->SetBranchAddress( "costheta2", &costheta2 );
  chainDT->SetBranchAddress( "costheta3", &costheta3 );
  chainDT->SetBranchAddress( "m12", &m12 );
  chainDT->SetBranchAddress( "m13", &m13 );
  chainDT->SetBranchAddress( "m23", &m23 );
 
  entriesMC = (Int_t)chainMC->GetEntries();
  entriesDT = (Int_t)chainDT->GetEntries();
 
  m12_low = chainDT->GetMinimum( "m12" );
  m13_low = chainDT->GetMinimum( "m13" );
  m23_low = chainDT->GetMinimum( "m23" );
 
  m12_up = chainDT->GetMaximum( "m12" );
  m13_up = chainDT->GetMaximum( "m13" );
  m23_up = chainDT->GetMaximum( "m23" );
 
  Int_t ny = calculateBins( entriesDT ); // bins in mass axisis (Y-axisis)
 
  MC1 = new TH2F( "myMC1", "", Ncos, -1.0, 1.0, ny, m23_low, m23_up ); // costheta1 vs. m23
  MC2 = new TH2F( "myMC2", "", Ncos, -1.0, 1.0, ny, m13_low, m13_up ); // costheta2 vs. m13
  MC3 = new TH2F( "myMC3", "", Ncos, -1.0, 1.0, ny, m12_low, m12_up ); // costheta3 vs. m12
 
  DT1 = new TH2F( "myDT1", "", Ncos, -1.0, 1.0, ny, m23_low, m23_up ); // costheta1 vs. m23
  DT2 = new TH2F( "myDT2", "", Ncos, -1.0, 1.0, ny, m13_low, m13_up ); // costheta2 vs. m13
  DT3 = new TH2F( "myDT3", "", Ncos, -1.0, 1.0, ny, m12_low, m12_up ); // costheta3 vs. m12
 
  WT1 = new TH2F( "myWT1", "", Ncos, -1.0, 1.0, ny, m23_low, m23_up ); // costheta1 vs. m23
  WT2 = new TH2F( "myWT2", "", Ncos, -1.0, 1.0, ny, m13_low, m13_up ); // costheta2 vs. m13
  WT3 = new TH2F( "myWT3", "", Ncos, -1.0, 1.0, ny, m12_low, m12_up ); // costheta3 vs. m12
 
  MC1->SetDirectory( 0 );
  MC2->SetDirectory( 0 );
  MC3->SetDirectory( 0 );
 
  DT1->SetDirectory( 0 );
  DT2->SetDirectory( 0 );
  DT3->SetDirectory( 0 );
 
  WT1->SetDirectory( 0 );
  WT2->SetDirectory( 0 );
  WT3->SetDirectory( 0 );
 
  // filling MC histogram
  for ( Int_t j = 0; j < entriesMC; j++ )
  {
    chainMC->GetEntry( j );
    MC1->Fill( costheta1, m23 );
    MC2->Fill( costheta2, m13 );
    MC3->Fill( costheta3, m12 );
  }
 
  // filling data histogram
  for ( Int_t j = 0; j < entriesDT; j++ )
  {
    chainDT->GetEntry( j );
    DT1->Fill( costheta1, m23 );
    DT2->Fill( costheta2, m13 );
    DT3->Fill( costheta3, m12 );
  }
  //--------debugging ---------
  //  showFrame(MC1);
  //  showFrame(MC2);
  //  showFrame(MC3);
 
  EvtHis2F::init( MC1, DT1, WT1 );
  WT1  = EvtHis2F::getHwt();
  max1 = EvtHis2F::getZmax();
 
  EvtHis2F::init( MC2, DT2, WT2 );
  WT2  = EvtHis2F::getHwt();
  max2 = EvtHis2F::getZmax();
 
  EvtHis2F::init( MC3, DT3, WT3 );
  WT3  = EvtHis2F::getHwt();
  max3 = EvtHis2F::getZmax();
 
  //------------
 
  /*
   show(WT1);
   std::cout<<"================================================="<<std::endl;
   show(WT2);
   std::cout<<"================================================="<<std::endl;
   show(WT3);
  */
}
 
int EvtNT3::calculateBins( int entries ) {
  int bins, ncell;
  ncell = 30; // at lease to require each cell to have 30 events
  bins  = entries / ncell / Ncos;
  if ( bins > 100 ) { bins = 100; }
  return bins;
}
 
bool EvtNT3::AR1( double costheta, double mass ) {
  bool accept = false;
  accept      = EvtHis2F::AR( costheta, mass, max1, WT1 );
  //--- debugging
  //  std::cout<<"Max_mass= "<<getZmax()<<std::endl;
 
  return accept;
}
 
bool EvtNT3::AR2( double costheta, double mass ) {
  bool accept = false;
  accept      = EvtHis2F::AR( costheta, mass, max2, WT2 );
 
  return accept;
}
 
bool EvtNT3::AR3( double costheta, double mass ) {
  bool accept = false;
  accept      = EvtHis2F::AR( costheta, mass, max3, WT3 );
 
  return accept;
}