#include 
#include 
#include 
#include "ntpana.h"
#include "gen_part.h"
#include "reco_part.h"
#include "cluster.h"
#include "pair.h"
#include "fourVector.h"
#include "ntputy.h"
#include "cxx_hbook.h"
#include "hrf_util.h"
#include "d0prop.h"
#include "runlist.h"
#include "header.h"
#include "dstar.h"
#include "dzero.h"
#include "lumi.h"
#include "bestDz.h"
#include "fisher.h"
#include "usefsh.h"

// Fisher analysis
const int nvar(3) ;
int Status(1);                         // Status: possible values 0, 1 , 2; indicates the step of Fisher analysis
Fisher fsh( nvar ) ;
float tag[nvar];

//#include "vub.h"
// shared variables
int nprtDs = 0 ;
int nprtDz = 0 ;
int nprtDcy= 0 ;
int nevt = 0   ; // event counter
int Bad  = 0   ; // rejected by runqualy
int OnPeak = 0 ; // on peak
int OffPeak = 0; // off peakint
int isData     ; // 0/1 (MC/data)
int ntpId      ; // ntuple address
int ntpRecSize ; // ntuple record size
int Skim=2     ; // input/output compressed ntuple
int nSkim = 0 ;
int Wgt = 1    ; // normalisation weight
int year;
int Sample ;
float Wlumi = 1    ; // normalisation weight
int hb(0),idh(0) ;

int* nRun;
int* nPair;
int* nTrack;
int* nClus;
int* nMC ;

char* chKind = new char [100] ;

const int ntag = 5 ;
const int nCones = 5  ; // cone cuts
const int Samples = 4 ; // mass/side good/wrong

const int TagN2 = 23 ;
const int TagN1 = 22 ;
const int TagKa = 13 ;
const int TagMu = 15 ;
const int TagEl = 14 ;
const int TagKm = 12 ;
const int TagKe = 11 ;

// ---------------------------------------------
// output units: one per each cone per four
//               samples times 2 (on/off peak)
// ---------------------------------------------

const int Units = Samples * nCones ;
ofstream OutOnPeak[Units] ;
ofstream OutOffPeak[Units] ;
ofstream OutMCSim[Units];

float KaMass = 0.494 ; // sara' giusta ?
float K0Mass = 0.498 ; // sara' giusta ?
float PiMass = 0.139 ;
float MuMass = 0.115 ;
float ElMass = 0.0005 ;
float cSpeed = 0.02997 ; // cm/ps

// tag variables
float pLepMin = 1.5  ;
float pLepMax = 2.5  ;
float pPstMin = 0.05 ;
float pPstMax = 0.17 ;
float VtxProb = 0.1  ;
float TagProb = 0.   ;
int   TagMult = 0 ;

const float LowMassBand = -2.0 ;
const float HigMassBand =  2.5 ;
const float LowSideBand =-10.0 ;
const float HigSideBand = -4.0 ;

const float cutEz = 0.3 ; // maximum error on Dz

// internal packages
void AnaIni() ;
void AnaEnd() ;
void OpenOutFiles();
void CloseOutFiles();
void WriteVtx( int , int , int , Pair* , int* , int* ) ;
void book_pair( int, int ) ;
void fill_pair( Pair*, int ,int* ) ;
void fill_by_kind( Pair*, int, int* ) ;
void fill_vertex( Pair*, int ,int* ) ;
void book_vertex( int, int ) ;
int EncodeDss( GenPart ) ;
// interface to Fortran Ntuple in-out package
extern "C" int ntpsteer_( int &, int &, int &, int & ) ;
extern "C" int noutend_( int & ) ;
extern "C" int noutwrt_dt_( int & ) ;
extern "C" int noutini_dt_( int & , int & ) ;
extern "C" void binout_mc_() ;
extern "C" void bininp_mc_() ;
extern "C" void binout_dt_() ;
extern "C" void bininp_dt_() ;
extern "C" void binend_() ;
extern "C" void binread_( int& ) ;
extern "C" void binini_( int &  ) ;

// main program //
int main(){

  // initialisation
  AnaIni() ;
  fshini( Status, &fsh ) ;

  // steering fortran from ntpana_ftn.o
  if( Skim == 2 )
    { binread_( isData ) ; }
  else
    {
      int Dum = 0;
      ntpsteer_( ntpId, ntpRecSize , isData, Dum ) ;
    }
  // end of analysis
  AnaEnd() ;
  fshend( Status, &fsh ) ;

  return 0 ;

}

void AnaEnd() {}


/*  ==================================================
    Initialisation: read title file to drive the run
    ================================================== */

void  AnaIni(){

  ifstream* InpTitles = new ifstream( "titles.data" ) ;

  // data or Monte Carlo ?
  cout << " enter 0 if Monte Carlo, else 1 " ;
  *InpTitles  >> isData ;  cout << endl ;

  // ntuple address
  cout << " enter NTP id number " ;
  *InpTitles  >> ntpId ;   cout << endl ;

  // get size of the input ntuple
  cout << " enter NTP Record Size " ;
  *InpTitles  >> ntpRecSize ;   cout << endl ;

  // event kind
  cout << " enter event kind " ; *InpTitles >> chKind ;
  cout << chKind << endl ;

  // year
  cout << " enter processing year " << endl;
  *InpTitles  >> year ;   cout << endl ;

  // luminosity weight
  float Lum(0) ;
  Sample = SampleKind( year,chKind,Lum,Wlumi ) ;
  hb = 1000*Sample ;
  delete InpTitles ;

  cout << "data Kind       " << chKind << " " << year << endl;
  cout << " Luminosity " << Lum << " Weight " << Wlumi << endl;

 cout << " *************************************** " << endl;
 cout << " OutPut Files For Analysis: " << endl << endl;
}

/*   ===============================================

     Analysis Skeleton. Called once per
     event by NTLOAD from Patch ntpana_main.car
     (also in ntpana_ftn.obj)

     =============================================== */


extern "C" int ntpana_( int *nRun, int *nPair, int *nTrack , int *nClus, int* nDz , int* nDst , int *nMC ) {

  if( *nTrack > 99 ) return 0 ;
  nevt++;
  if( nevt%1000 == 0 ){cout << endl
			    << " Evt. " << nevt
			    << " Skim " << nSkim < -2
     Side Band -8 +-+ -3 (id+50)
     otherwise skip
  */

  if( nuMass < -2 ) return 0;

  fourVector B04V = thePair.B04V();
  fourVector Pst4V = thePst.Rf4V();
  fourVector Dst4V = thePair.Dst4V();
  fourVector lTag4V = thePair.Lep4V();

  // loop over tracks;
  // get good charge electrons and kaons

  float QPst = thePst.Charge() ;

  int nLep(0),LepList[10] ;
  int nKch(0),KchList[10] ;

  for( int i=0 ; i<*nTrack ; i++ ) {

    if( i == thePair.LepId() ) continue ;
    if( i == thePair.PstId() ) continue ;
    RecoPart aPart = RecoPart( nTrack,i ) ;

    //    if( aPart.TagEl() && aPart.Charge() == QPst )
    if( aPart.TagEl() )
      {
	LepList[nLep]=i;
	nLep++;
      }
    else if (aPart.TagKa() && aPart.Charge() == -QPst )
      {
	KchList[nKch]=i;
	nKch++;
      }
  }

  int DecayMode(-1) ;
  int GenDzId(-1),GenDzLep(-1), GenDzHad(-1), GenDzNeu(-1) ;
  int GenDstId(-1);

  float QsqLep(-1),PgenLep(-1) ;
  float QsqHad(-1),PgenHad(-1) ;
  fourVector DzGen4V,DstGen4V ;


  if( thePair.IsSignal() )
    {
      GenDzId = thePair.D0GenAdd() ;
      GenDstId = thePair.DstGenAdd() ;
      if( GenDzId > 0 )
	{
	  theDzGen = GenPart( nMC, GenDzId ) ;
	  DzGen4V = theDzGen.Rf4V() ;
	  DecayMode = DztoElX( theDzGen, nMC,
			       GenDzLep,GenDzHad,GenDzNeu );
	}

      Pawc.hfill( 91,DecayMode,0.,Wlumi ) ;
      if( DecayMode > 0 ) {

	int ListDz[10];
	int mult = theDzGen.ListDau( ListDz );

	//	if( DecayMode == 1 ) theDzGen.DumpList( mult,ListDz ) ;

	int jd = 100*DecayMode ;

	Pawc.hfill( jd+1, mult , 0., Wlumi ) ;

	GenPart theHad = GenPart( nMC,GenDzHad ) ;
	GenPart theEle = GenPart( nMC,GenDzLep ) ;
	GenPart theNeu = GenPart( nMC,GenDzNeu ) ;

	fourVector Had4V = theDzGen.Rf4V()-theHad.Rf4V() ;
	fourVector Lep4V = theEle.Rf4V()+theNeu.Rf4V() ;
	QsqLep = Lep4V.mSqr() ;
	QsqHad = Had4V.mSqr() ;
	PgenLep = theEle.PmagRf() ;
        PgenHad = theHad.PmagRf() ;
	Pawc.hfill( jd+2, QsqLep , QsqHad , 1. ) ;
	if( mult == 3 ) Pawc.hfill( jd+3, QsqLep , QsqHad , 1. ) ;
	Pawc.hfill( jd+4, PgenLep , 0.,1. ) ;
	Pawc.hfill( jd+5, PgenHad , 0.,1. ) ;

      }
    }

  // fill measured quantities

  if( nLep == 0 ) return 0 ;
  if( nKch == 0 ) return 0 ;

  Pawc.hfill( idh+902, nuMass, 0., 1.) ;
  if( Skim==1 ) {
    nSkim++;
    if( isData )
      { binout_dt_() ; }
    else
      { binout_mc_();}
  }

  int evmul(0);
  for( int iLep = 0 ; iLep 0.4) (varing_M= 0.3999999);
	      if( Status == 2 ) Pawc.hfill( 666012, varing_M, 0, Wlumi ) ;

	    }
	//fine montecarlo
      if( !isData )
	{

	  //.  float tag[nvar];
	  tag[0] = neu4V.mSqr();
	  tag[1] = DM ;
	  tag[2] = MD ;
	  float varing_S(0) ;
	  float varing_B(0) ;

	  if( DecayMode == 1 ) // segnale
	    {
	      varing_S= fshana( Status, 1 , &fsh , nvar , tag ) ;
	      if (varing_S< -0.2) (varing_S= -0.1999999);
	      if (varing_S> 0.4) (varing_S= 0.3999999);
	      if( Status == 2 ) Pawc.hfill( 13001, varing_S, 0, Wlumi ) ;
	    }
	  else if( DecayMode == 0 || DecayMode == 6 ) //fondo
	    {
	      varing_B= fshana (Status , 0 , &fsh , nvar , tag );
	      if (varing_B< -0.2) (varing_B= -0.1999999);
	      if (varing_B> 0.4) (varing_B= 0.3999999);
	      if( Status == 2 ) Pawc.hfill( 14001, varing_B, 0, Wlumi) ;
	    }
	}
      evmul++ ;
    }
  }
  Pawc.hfill( idh+10, evmul , 0. , Wlumi ) ;

  return 0 ;
} // end ntpana



extern "C" void anaboo_() {

  if( Skim == 1 ) {
    int iread(0);
    if( isData )
      { binini_(iread) ; }
    else
      { binini_(iread) ; }
  }

  // control histos
  Pawc.hbook1( 2, " CME  ",100,9.,11.,0.) ;
  Pawc.hbook1( 3, " [bg] ",50, 0.5,0.6,0.)  ;
  Pawc.hbook1( 4, " flavour ",20,510.,530.,0.) ;
  Pawc.hbook1( 8, " Lepton Tag ",50,0.,50.,0.) ;
  Pawc.hbook1( 9, " Lepton Tag ",50,0.,50.,0.) ;
  Pawc.hbook1( 91," D^0! decay mode ", 11,0.,11.,0. ) ;

  Pawc.hbook1( hb+901 , " M?[n,1]^2! ", 65,-10.,3.,0. ) ;
  Pawc.hbook1( hb+902 , " M?[n,1]^2! ", 65,-10.,3.,0. ) ;

  if( !isData )
    {
      Pawc.hbook1( hb+2901 , " M?[n,1]^2! ", 65,-10.,3.,0. ) ;
      Pawc.hbook1( hb+2902 , " M?[n,1]^2! ", 65,-10.,3.,0. ) ;
    }

  for( int js = -1 ; js < 2 ; js+=2 ) {
    for( int ks = 0 ; ks <100 ; ks+=50 ) {

      int kb = hb+ks ;
      Pawc.hbook1( js*(kb+10 ), " Mult " , 5,0.,5.,0. )  ;
      Pawc.hbook2( js*(kb+1 ), "P vs [Q] K ", 50,0.,1., 50,0.,2.5, 0. ) ;
      Pawc.hbook2( js*(kb+2 ), "P vs [Q] e ", 50,0.,1., 50,0.,2.5, 0. ) ;
      Pawc.hbook2( js*(kb+3 ), "[Q]?K! vs [Q]?e! ", 50,0.,1., 50,0.,1.,0. ) ;
      Pawc.hbook1( js*(kb+4 ), "[n] m ^2! ", 65,-10.,3.,0.);
      Pawc.hbook2( js*(kb+5 ), "Q2 vs [n] m ^2! ", 45,-6.,3.,10,0.,3.,0.);
      Pawc.hbook2( js*(kb+7 ), "[D]M vs [n]M^2! ",30,-4.5,1.5,60,0.130,0.190,0. ) ;
      Pawc.hbook2( js*(kb+8 ), "[D]M cns vs [n]M^2! ",30,-4.5,1.5,30,0.140,0.170,0. ) ;
      Pawc.hbook2( js*(kb+6 ), "Mass vs [n]M^2! ",30,-4.5,1.5,40,0.,2.,0. ) ;
      Pawc.hbprof( js*(kb+9 ), "[D]M vs [n]M^2! ",60,-4.5,1.5,0.130,0.190," " ) ;
      if( !isData ) {
	Pawc.hbook2( js*(kb+2001 ),"P vs [Q] K ", 50,0.,1., 50,0.,2.5, 0. ) ;
	Pawc.hbook2( js*(kb+2002 ),"P vs [Q] e ", 50,0.,1., 50,0.,2.5, 0. ) ;
	Pawc.hbook2( js*(kb+2003 ),"[Q]?K! vs [Q]?e! ", 50,0.,1., 50,0.,1.,0. ) ;
	Pawc.hbook1( js*(kb+2004 ),"[n] m ^2! ", 65,-10.,3.,0.);
	Pawc.hbook2( js*(kb+2005 ),"Q2 vs [n] m ^2! ", 45,-6.,3.,10,0.,3.,0.);
	Pawc.hbook2( js*(kb+2007 ),"[D]M vs [n]M^2! ",30,-4.5,1.5,60,0.130,0.190,0. ) ;
	Pawc.hbook2( js*(kb+2008 ),"[D]M cns vs [n]M^2! ",30,-4.5,1.5,30,0.140,0.170,0.);
	Pawc.hbook2( js*(kb+2006 ),"Mass vs [n]M^2! ",30,-4.5,1.5,40,0.,2.,0. ) ;
	Pawc.hbprof( js*(kb+2009 ), "[D]M vs [n]M^2! ",60,-4.5,1.5,0.130,0.190," " ) ;
      }
    }
  }

  if( isData ) return ;

  for( int id = 100 ; id<700 ; id += 50 ) {
    Pawc.hbook1( id+1 , " mult " , 10, 0., 10.,0. );
    Pawc.hbook2( id+2 , "Qhad vs Qlep " , 60,0.,3., 60,0.,3. , 0.) ;
    Pawc.hbook2( id+3 , "Qhad vs Qlep " , 60,0.,3., 60,0.,3. , 0.) ;
    Pawc.hbook1( id+4 , "P e" , 50, 0.,2.5 ,0. ) ;
    Pawc.hbook1( id+5 , "P k" , 50, 0.,2.5 ,0. ) ;
    Pawc.hbook1( id+6 , "P e sel" , 50, 0.,2.5 ,0. );
    Pawc.hbook1( id+7 , "P k sel" , 50, 0.,2.5 ,0. );
    Pawc.hbook1( id+8 , "[D]Q" , 50, -2.5,2.5 ,0. );
    Pawc.hbook2( id+9 , "[D]Q vs Q?gen!" ,30,0.,3.,25,-2.5,2.5,0.);
    Pawc.hbook2( id+10 , "[D]Q vs Q?rec!" ,30,0.,3.,25,-2.5,2.5,0.);
    Pawc.hbprof( id+11 , "[D]Q vs Q?rec!" ,30,0.,3.,-2.5,2.5," ");
    Pawc.hbook1( 666012, "Confronto Montecalro", 100, -0.2, 0.4, 0. );

    Pawc.hbook1( id+20 , "[D]E" , 50, -1.5,1.5 ,0. );
    Pawc.hbook2( id+21 , "[D]E vs E" ,50,1.5,3.5,25,-1.0,1.0,0. );
    Pawc.hbprof( id+22 , "[D]E vs E" ,50,1.5,3.5,-1.0,1.0," " );
    Pawc.hbook1( id+23 , " cos([Q])" , 50, 0.,1. ,0. );
    Pawc.hbook2( id+24 , "[D]Q vs [D]E" ,25,-1.0,1.0,25,-2.5,2.5,0.);
    Pawc.hbprof( id+25 , "[D]Q vs [D]E" ,25,-1.0,1.0,-2.5,2.5," ");
    Pawc.hbook2( id+26 , "[D]Q vs [D]E" ,30,0.7,1.0,25,-2.5,2.5,0.);

    Pawc.hbook2( hb+id+1, "P vs [Q] K ", 50,0.,1., 50,0.,2.5, 0. ) ;
    Pawc.hbook2( hb+id+2, "P vs [Q] e ", 50,0.,1., 50,0.,2.5, 0. ) ;
    Pawc.hbook2( hb+id+3, "[Q]?K! vs [Q]?e! ", 50,0.,1., 50,0.,1.,0. ) ;
    Pawc.hbook1( hb+id+4, "[n] m ^2! ", 65,-10.,3.,0.);
    Pawc.hbook2( hb+id+5, "Q2 vs [n] m ^2! ", 45,-6.,3.,10,0.,3.,0.);
    Pawc.hbook2( hb+id+7, "[D]M vs [n]M^2! ",30,-4.5,1.5,60,0.130,0.190,0. ) ;
    Pawc.hbook2( hb+id+8, "[D]M cns vs [n]M^2! ",30,-4.5,1.5,30,0.140,0.170,0. ) ;
    Pawc.hbook2( hb+id+6, "Mass vs [n]M^2! ",30,-4.5,1.5,40,0.,2.,0. ) ;
    Pawc.hbprof( hb+id+9, "[D]M vs [n]M^2! ",60,-4.5,1.5,0.130,0.190," " ) ;

    Pawc.hbook1( 13001, "Signal", 100, -0.2, 0.4,0.) ;
    Pawc.hbook1( 14001, "BacK", 100, -0.2, 0.4, 0. );
    Pawc.hbook1( 666012, "Confronto Montecalro", 100, -0.2, 0.4, 0. );



  }
  Pawc.hbarx(0) ;

  return;
}



extern "C" void anaend_() {

  /*  if( Skim == 1 ) {
    int ntpout = 111;
    noutend_(ntpout);
    }*/

  /*
    HRFile hrF ( "ntpana.hist" , New ) ;
    int icycle = 0 ;
    hrF.hrout(  Pawc,0,icycle," " ) ;
  */

  binend_() ;
  return;

}

    Source: geocities.com/it/realquasars/stage

               ( geocities.com/it/realquasars)                   ( geocities.com/it)