/////////////////////////////////////////////////////////
//   This class has been automatically generated 
//     (Mon Sep 24 10:55:41 2001 by ROOT version3.01/06)
//   from TTree tuple/Selex reconstruction data
//   found on file: /export/data/ayan/passII_root_data/oc0_vee_1.root
//////////////////////////////////////////////////////////


#ifndef s_h
#define s_h

#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>

class s {
   public :
   TTree          *fChain;   //pointer to the analyzed TTree or TChain
   Int_t           fCurrent; //current Tree number in a TChain
//Declaration of leaves types
   Float_t         RN_run;
   Float_t         RN_event;
   Float_t         RN_ID;
   Float_t         RN_x0b;
   Float_t         RN_y0b;
   Float_t         RN_Pxb;
   Float_t         RN_Pyb;
   Float_t         RN_Pzb;
   Float_t         RN_Pidb;
   Float_t         RN_Xp;
   Float_t         RN_Yp;
   Float_t         RN_Zp;
   Float_t         RN_Sxp;
   Float_t         RN_Syp;
   Float_t         RN_Szp;
   Float_t         RN_Pchi;
   Float_t         RN_Pnf;
   Float_t         RN_Pnt;
   Float_t         RN_mass;
   Float_t         RN_Pxs;
   Float_t         RN_Pys;
   Float_t         RN_Pzs;
   Float_t         RN_Es;
   Float_t         RN_Xs;
   Float_t         RN_Ys;
   Float_t         RN_Zs;
   Float_t         RN_Sxs;
   Float_t         RN_Sys;
   Float_t         RN_Szs;
   Float_t         RN_Schi;
   Float_t         RN_Isol;
   Float_t         RN_mass2;
   Float_t         RN_Pxs2;
   Float_t         RN_Pys2;
   Float_t         RN_Pzs2;
   Float_t         RN_Es2;
   Float_t         RN_Xs2;
   Float_t         RN_Ys2;
   Float_t         RN_Zs2;
   Float_t         RN_Sxs2;
   Float_t         RN_Sys2;
   Float_t         RN_Szs2;
   Float_t         RN_Schi2;
   Float_t         RN_Isol2;
   Float_t         RN_tgt;
   Float_t         RN_nt;
   Float_t         RN_ns;
   Float_t         RN_Xf;
   Float_t         RN_Pt;
   Float_t         RN_L;
   Float_t         RN_S;
   Float_t         RN_LoS;
   Float_t         RN_pvtx;
   Float_t         RN_smin;
   Float_t         RN_scut;
   Float_t         RN_X0t1;
   Float_t         RN_Y0t1;
   Float_t         RN_Pxt1;
   Float_t         RN_Pyt1;
   Float_t         RN_Pzt1;
   Float_t         RN_chit1;
   Float_t         RN_masst1;
   Float_t         RN_pidt1;
   Float_t         RN_pmtrkt1;
   Float_t         RN_pid1t1;
   Float_t         RN_pid2t1;
   Float_t         RN_mvert_t1;
   Float_t         RN_X0t2;
   Float_t         RN_Y0t2;
   Float_t         RN_Pxt2;
   Float_t         RN_Pyt2;
   Float_t         RN_Pzt2;
   Float_t         RN_chit2;
   Float_t         RN_masst2;
   Float_t         RN_pidt2;
   Float_t         RN_pmtrkt2;
   Float_t         RN_pid1t2;
   Float_t         RN_pid2t2;
   Float_t         RN_mvert_t2;
   Float_t         RN_X0t3;
   Float_t         RN_Y0t3;
   Float_t         RN_Pxt3;
   Float_t         RN_Pyt3;
   Float_t         RN_Pzt3;
   Float_t         RN_chit3;
   Float_t         RN_masst3;
   Float_t         RN_pidt3;
   Float_t         RN_pmtrkt3;
   Float_t         RN_pid1t3;
   Float_t         RN_pid2t3;
   Float_t         RN_mvert_t3;
   Float_t         RN_X0t4;
   Float_t         RN_Y0t4;
   Float_t         RN_Pxt4;
   Float_t         RN_Pyt4;
   Float_t         RN_Pzt4;
   Float_t         RN_chit4;
   Float_t         RN_masst4;
   Float_t         RN_pidt4;
   Float_t         RN_pmtrkt4;
   Float_t         RN_pid1t4;
   Float_t         RN_pid2t4;
   Float_t         RN_mvert_t4;
   Float_t         RN_X0t5;
   Float_t         RN_Y0t5;
   Float_t         RN_Pxt5;
   Float_t         RN_Pyt5;
   Float_t         RN_Pzt5;
   Float_t         RN_chit5;
   Float_t         RN_masst5;
   Float_t         RN_pidt5;
   Float_t         RN_pmtrkt5;
   Float_t         RN_pid1t5;
   Float_t         RN_pid2t5;
   Float_t         RN_mvert_t5;
   Float_t         RN_X0t6;
   Float_t         RN_Y0t6;
   Float_t         RN_Pxt6;
   Float_t         RN_Pyt6;
   Float_t         RN_Pzt6;
   Float_t         RN_chit6;
   Float_t         RN_masst6;
   Float_t         RN_pidt6;
   Float_t         RN_pmtrkt6;
   Float_t         RN_pid1t6;
   Float_t         RN_pid2t6;
   Float_t         RN_mvert_t6;
   Float_t         RN_X0t7;
   Float_t         RN_Y0t7;
   Float_t         RN_Pxt7;
   Float_t         RN_Pyt7;
   Float_t         RN_Pzt7;
   Float_t         RN_chit7;
   Float_t         RN_masst7;
   Float_t         RN_pidt7;
   Float_t         RN_pmtrkt7;
   Float_t         RN_pid1t7;
   Float_t         RN_pid2t7;
   Float_t         RN_mvert_t7;
   Float_t         RN_X0t8;
   Float_t         RN_Y0t8;
   Float_t         RN_Pxt8;
   Float_t         RN_Pyt8;
   Float_t         RN_Pzt8;
   Float_t         RN_chit8;
   Float_t         RN_masst8;
   Float_t         RN_pidt8;
   Float_t         RN_pmtrkt8;
   Float_t         RN_pid1t8;
   Float_t         RN_pid2t8;
   Float_t         RN_mvert_t8;
   Float_t         RN_X0t9;
   Float_t         RN_Y0t9;
   Float_t         RN_Pxt9;
   Float_t         RN_Pyt9;
   Float_t         RN_Pzt9;
   Float_t         RN_chit9;
   Float_t         RN_masst9;
   Float_t         RN_pidt9;
   Float_t         RN_pmtrkt9;
   Float_t         RN_pid1t9;
   Float_t         RN_pid2t9;
   Float_t         RN_mvert_t9;
   Float_t         RN_chi_rv02;
   Float_t         RN_dca_rv02;
   Float_t         RN_tet_rv02;
   Float_t         RN_pd_ph;
   Float_t         RN_chi_rv01;
   Float_t         RN_dca_rv01;
   Float_t         RN_pidtrk_rv01;
   Float_t         RN_pd_ph_b;
   Float_t         RN_chi_rv0b;
   Float_t         RN_dca_rv0b;
   Float_t         RN_z_rv0b;
   Float_t         RN_vid_rv01b;
   Float_t         RN_vid_rv02b;
   Float_t         RN_status_rv0;

//List of branches
   TBranch        *b_RN;

   char decayid[3];
   int idecay;
   char ps_filename[80];
   char pavetext[256];
   char *cid [9];
   Float_t kaonprob[4];
   Float_t pionprob[4];
   
   s(TTree *tree=0);
   ~s();
   Int_t  Cut(Int_t entry);
   Int_t  GetEntry(Int_t entry);
   Int_t  LoadTree(Int_t entry);
   
   Int_t TgtVol(Int_t entry);        // new 6/19/01 L. Dauwe
   Int_t PzBeam(Int_t entry);        // new 6/15/01 L. Dauwe
   Int_t PartidBeam(Int_t entry);    // new 6/22/01 L. Dauwe
   Int_t Daughter_Is_A_Pion(Int_t entry);
   Int_t Daughter_Is_A_Kaon(Int_t entry);
   Int_t Mvert_Part(Int_t entry, Int_t ver_indx, Float_t Momentum, const char *part_name);
   Int_t DaughterKaonProb(Int_t entry);
   Int_t DaughterPionProb(Int_t entry);
   Int_t MinOfTrk2Trk3(Int_t entry, float mom_thres );
   Int_t MinOfTrk2(Int_t entry, float mom_thres);
   void kp(Int_t entry);
   Int_t Event_Reject(Int_t entry);
   Int_t TrkChi2(Int_t entry, Int_t ver_indx); 
   Double_t Trk_p(Int_t entry, Int_t trk); // calculate track momentum
   
   Double_t Oi_p_asym(Int_t entry); // 850 Ome- pi+ momentum asymmetry
   
   Double_t oiiipt(Int_t entry);
   Double_t oipt(Int_t entry);
   Double_t sum_pt2(Int_t entry);
   Double_t sum_pt(Int_t entry);
   Double_t pts(Int_t entry);

   Float_t smin(Int_t entry, Int_t position); 
   void Lambda0(Int_t entry) ;
   Int_t Angle_Vector(Int_t entry);
   Double_t ProperTime(Int_t entry,Int_t loversigma);
   Double_t Cos_Theta(Int_t entry);
   Int_t kk827(Int_t entry) ;
   Int_t pid828(Int_t entry);
   Int_t pid829(Int_t entry);
   Int_t pid850(Int_t entry);
   Int_t BeamId(Int_t entry);

   Int_t VidPionProb(Int_t entry) ;
   Int_t Vid_Is_APion(Int_t entry); 
   Int_t VidKaonProb(Int_t entry) ;
   Int_t Vid_Is_AKaon(Int_t entry); 
 
   Double_t ReCalc4(Int_t entry);
   Double_t ReCalc2(Int_t entry);
   Double_t NewMass(Int_t entry);
   Double_t Drmass(Int_t entry);

   Double_t DaughterMom(Int_t entry);

   void Slice2DHisto(TH2F *hp);
 
   void Order(int *p,int *q);
   void Bubble(int *a,int n);
   void hex2bin(Int_t n);
   Int_t FindParticle(Int_t entry, const char *part_name);

   void GenerateEmbedFile(Int_t entry);
   void DumpEvent(Int_t entry);
   void  NoteEvent(Int_t entry);
   void Init(TTree *tree);
   void Loop();
   void Show(Int_t entry = -1);
   void Track_Print(Int_t entry);
   void Histo_Normalizer();
   void ProperTimePlot();
   void OmeZs2ProTime();
   Bool_t Notify();
   
};

#endif

#ifdef s_cxx
s::s(TTree *tree)
{
// if parameter tree is not specified (or zero), connect the file
// used to generate this class and read the Tree.
   if (tree == 0) {


#ifdef SINGLE_TREE
      // The following code should be used if you want this class to access
      // a single tree instead of a chain
      TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("Memory Directory");
      if (!f) {
         f = new TFile("Memory Directory");
         f->cd("Rint:/");
      }
      tree = (TTree*)gDirectory->Get("tuple");

#else // SINGLE_TREE

      // The following code should be used if you want this class to access a chain
      // of trees.
          TChain * chain = new TChain("tuple","");
          //chain->Add("/home/ayan/passIIdata/p2b01.root");
          //chain->Add("/home/ayan/passIIdata/p2b02.root");
          //chain->Add("/home/ayan/passIIdata/p2c01.root");
          //chain->Add("/home/ayan/passIIdata/p2d01.root");
          //chain->Add("/home/ayan/passIIdata/p2e01.root");
          //chain->Add("/home/ayan/passIIdata/p2f01.root");
          //chain->Add("/home/ayan/passIIdata/p2g01.root");
          //chain->Add("/home/ayan/passIIdata/p2h01.root");
          //chain->Add("/home/ayan/passIIdata/p2p01.root");
          //chain->Add("/home/ayan/passIIdata/p2x01.root");
          //chain->Add("/home/ayan/passIIdata/p2z02.root");
          //chain->Add("/home/ayan/passIIdata/p2u01.root");//unfiltered
          //chain->Add("/home/ayan/passIIdata/p2z01.root");//same runs as p2z02
	  
	  //Monte Carlo
	  chain->Add("/home/ayan/Monte_Carlo/xf0/500K/oc0_oiii-xf0.root1");
	  chain->Add("/home/ayan/Monte_Carlo/xf0/500K/oc0_oi-xf0.root1");
	  //chain->Add("/home/ayan/Monte_Carlo/xf0/500K/oc0_xkii-xf0.root");
	  //chain->Add("/home/ayan/Monte_Carlo/xf3/500K/oc0_oiii_l0k-_xf3.root1");
	  //chain->Add("/home/ayan/Monte_Carlo/xf3/500K/oc0_oiii_l0k-_xf3.root2");
	  //chain->Add("/home/ayan/Monte_Carlo/xf3/500K/oc0_oi_l0k-_xf3.root1");
	  //chain->Add("/home/ayan/Monte_Carlo/xf3/500K/oc0_oi_l0k-_xf3.root2");
	  
	  tree = chain;
#endif // SINGLE_TREE


   }
   Init(tree);
}

s::~s()
{
   if (!fChain) return;
   delete fChain->GetCurrentFile();
}

Int_t s::GetEntry(Int_t entry)
{
// Read contents of entry.
   if (!fChain) return 0;
   return fChain->GetEntry(entry);
}
Int_t s::LoadTree(Int_t entry)
{
// Set the environment to read one entry
   if (!fChain) return -5;
   Int_t centry = fChain->LoadTree(entry);
   if (centry < 0) return centry;
   if (fChain->IsA() != TChain::Class()) return centry;
   TChain *chain = (TChain*)fChain;
   if (chain->GetTreeNumber() != fCurrent) {
      fCurrent = chain->GetTreeNumber();
      Notify();
   }
   return centry;
}

void s::Init(TTree *tree)
{
//   Set branch addresses
   if (tree == 0) return;
   fChain    = tree;
   fCurrent = -1;
   fChain->SetMakeClass(1);

   fChain->SetBranchAddress("RN",&RN_run);
   Notify();
}

Bool_t s::Notify()
{
//   called when loading a new file
//   get branch pointers
   b_RN = fChain->GetBranch("RN");
   return kTRUE;
}

void s::Show(Int_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
   if (!fChain) return;
   fChain->Show(entry);
}


Int_t s::BeamId(Int_t entry)
{
  Int_t sigma; Int_t pro; Int_t pion;
  
  sigma  = (int(RN_Pidb)%1000)/100;
  //omega  = (int(RN_Pidb)%10)/1;
  pro = (int(RN_Pidb)%10000)/1000;
  pion   = (int(RN_Pidb)%1000000)/100000;

  //cout<<sigma<<" "<<pro<<" "<<pion<<endl;
  
  if (sigma  > pion)  return 21; // Sigma- Beam
  if (pion   > sigma) return 9;  // Pion- Beam
  else  return 14;               // proton Beam
}

Float_t s::smin(Int_t entry, Int_t position) { 
  // smin is the distance from the primary vertex to an extension of a track
  // participating in a secondary vertex.
  // example: smin[0] is smin for the track with the largest smin value.
  // example: smin[1] is smin for the track with the second largest smin value.
Float_t s[4], smin[4];

  s[0] = sqrt( pow(RN_X0t1 + (RN_Pxt1/RN_Pzt1)*(RN_Zs - RN_L) - RN_Xp,2) + 
	       pow(RN_Y0t1 + (RN_Pyt1/RN_Pzt1)*(RN_Zs - RN_L) - RN_Yp,2) );
  s[1] = sqrt( pow(RN_X0t2 + (RN_Pxt2/RN_Pzt1)*(RN_Zs - RN_L) - RN_Xp,2) + 
	       pow(RN_Y0t2 + (RN_Pyt2/RN_Pzt1)*(RN_Zs - RN_L) - RN_Yp,2) );
  s[2] = sqrt( pow(RN_X0t3 + (RN_Pxt3/RN_Pzt1)*(RN_Zs - RN_L) - RN_Xp,2) + 
	       pow(RN_Y0t3 + (RN_Pyt3/RN_Pzt1)*(RN_Zs - RN_L) - RN_Yp,2) );
  s[3] = sqrt( pow(RN_X0t4 + (RN_Pxt4/RN_Pzt1)*(RN_Zs - RN_L) - RN_Xp,2) + 
	       pow(RN_Y0t4 + (RN_Pyt4/RN_Pzt1)*(RN_Zs - RN_L) - RN_Yp,2) );
 if (entry==502816) cout<<s[0]<<s[1]<<s[2]<<s[3]<<"jentry"<<entry <<"smin[position]"<<smin[position]<<endl;
  for (Int_t i=0; i<=3; i++) {
    smin[i] = 0;
    Int_t k = 0;
    for (Int_t j=0; j<=3; j++) {
      if (s[j]>0 && s[j]>smin[i]) {
	k = j;
	smin[i] = s[j];
      }
    }
    if (k >= 0) {
      s[k] = -s[k];
    }
  }
 
 return smin[position];
}

Int_t s::Daughter_Is_A_Pion(Int_t entry) {  
 
Float_t DaughterPionProb=0.; // V0 Kink Pion Probability
Float_t DaughterKaonProb=0.; // V0 Kink Kaon Probability
Float_t DaughterProtonProb=0.; // V0 Kink Proton Probability

DaughterPionProb   = ( int(RN_pidtrk_rv01) % 100 )  ;
DaughterKaonProb   = ( int(RN_pidtrk_rv01)/100 % 100 ) ;
DaughterProtonProb = ( int(RN_pidtrk_rv01)/10000 % 100 ) ;

if (DaughterPionProb > DaughterProtonProb && DaughterPionProb > DaughterKaonProb) 
//if (DaughterPionProb > DaughterKaonProb) 
  {
    return 1 ; // Daughter IS a Pion
  }
else 
  {
    return -1;  // Daughter is NOT a Pion
  }
}

Int_t s::Daughter_Is_A_Kaon(Int_t entry) {  
  /* V0 Kink Kaon Probability*/
Float_t DaughterPionProb=0.;
Float_t DaughterKaonProb=0.;
Float_t DaughterProtonProb=0.;
 
DaughterPionProb   = ( int(RN_pidtrk_rv01) % 100 )  ;
//DaughterKaonProb   = ( int(RN_pidtrk_rv01)/100 % 100 ) ;
//DaughterProtonProb = ( int(RN_pidtrk_rv01)/10000 % 100 ) ;
//swap to correct it
DaughterKaonProb   = ( int(RN_pidtrk_rv01)/100 % 100 ) ;
DaughterProtonProb = ( int(RN_pidtrk_rv01)/10000 % 100 ) ;
 
if (DaughterKaonProb >DaughterProtonProb && DaughterKaonProb > DaughterPionProb) 
//if (DaughterKaonProb > DaughterPionProb) 
   {
     return 1 ; // Daughter IS a Pion
   }
 else 
   {
     return -1;  // Daughter is NOT a Pion
   }
}
Int_t s::DaughterKaonProb(Int_t entry) {  
  /* V0 Kink Kaon Probability*/
Float_t DaughterPionProb=0.;
Float_t DaughterKaonProb=0.;
Float_t DaughterProtonProb=0.;
DaughterPionProb   = ( int(RN_pidtrk_rv01) % 100 )  ;
DaughterKaonProb   = ( int(RN_pidtrk_rv01/100) % 100 ) ;
DaughterProtonProb = ( int(RN_pidtrk_rv01/10000) % 100 ) ;
float sig = ( int(RN_pidtrk_rv01/1000000) % 100 ) ;
Float_t ricid=  ( int(RN_pidtrk_rv01/100000000) % 10 )-1 ;

//cout<<ricid<<" "<<(RN_status_rv0 & ome_s_v0) <<" "<<int(RN_pidtrk_rv01)<<" "<<DaughterPionProb<<" "<<DaughterKaonProb<<" "<<DaughterProtonProb<<endl;
 if (DaughterKaonProb > 70) 
//if (DaughterKaonProb !=0) gives the same resulty as above
   {
     return 1 ; // Daughter IS a Kaon
   }
 else 
   {
     return -1;  // Daughter is NOT a Kaon
   }
}


Int_t s::DaughterPionProb(Int_t entry) {  
  /* V0 Kink Pion Probability*/
Float_t DaughterPionProb=0.;
Float_t DaughterKaonProb=0.;
Float_t DaughterProtonProb=0.;
DaughterPionProb   = ( int(RN_pidtrk_rv01) % 100 )  ;
DaughterKaonProb   = ( int(RN_pidtrk_rv01)/100 % 100 ) ;
DaughterProtonProb = ( int(RN_pidtrk_rv01)/10000 % 100 ) ;


 
 if (DaughterPionProb > 70) 
//if (DaughterPionProb !=0)  gives the same resulty as above
   {
     return 1 ; // Daughter IS a Pion
   }
 else 
   {
     return -1;  // Daughter is NOT a Pion
   }
}

Int_t s::VidPionProb(Int_t entry) {  
  /* V0 Kink Pion Probability vid_rv01*/
  double tkpiprob = (int(RN_vid_rv01b)%1000000-int(RN_vid_rv01b)%100000) / 1000000.;
  double tkkaprob = (int(RN_vid_rv01b)%100000 -int(RN_vid_rv01b)%10000)  / 100000.;
  double tkprprob = (int(RN_vid_rv01b)%10000  -int(RN_vid_rv01b)%1000)   / 10000.;
  if ( tkpiprob> .5) 
    {
      return 1 ; // Daughter IS a Pion
    }
  else 
    {
      return -1;  // Daughter is NOT a Pion
    }
}
Int_t s::Vid_Is_APion(Int_t entry) {  
  /* V0 Kink Pion Probability vid_rv01*/
  double tkpiprob = (int(RN_vid_rv01b)%1000000-int(RN_vid_rv01b)%100000) / 1000000.;
  double tkkaprob = (int(RN_vid_rv01b)%100000 -int(RN_vid_rv01b)%10000)  / 100000.;
  double tkprprob = (int(RN_vid_rv01b)%10000  -int(RN_vid_rv01b)%1000)   / 10000.;
  
  if ( tkpiprob>tkprprob && tkpiprob >tkkaprob ) return 1; 
  else return -1;

}

Int_t s::VidKaonProb(Int_t entry) {  
  /* V0 Kink Pion Probability vid_rv01*/
  double tkpiprob = (int(RN_vid_rv01b)%1000000-int(RN_vid_rv01b)%100000) / 1000000.;
  double tkkaprob = (int(RN_vid_rv01b)%100000 -int(RN_vid_rv01b)%10000)  / 100000.;
  double tkprprob = (int(RN_vid_rv01b)%10000  -int(RN_vid_rv01b)%1000)   / 10000.;
  if ( tkkaprob> .7) 
    {
      return 1 ; // Daughter IS a Pion
    }
  else 
    {
      return -1;  // Daughter is NOT a Pion
    }
}
Int_t s::Vid_Is_AKaon(Int_t entry) {  
  /* V0 Kink Pion Probability vid_rv01*/
  double tkpiprob = (int(RN_vid_rv01b)%1000000-int(RN_vid_rv01b)%100000) / 1000000.;
  double tkkaprob = (int(RN_vid_rv01b)%100000 -int(RN_vid_rv01b)%10000)  / 100000.;
  double tkprprob = (int(RN_vid_rv01b)%10000  -int(RN_vid_rv01b)%1000)   / 10000.;
  
  if ( tkkaprob>tkprprob && tkkaprob >tkpiprob ) return 1; 
  else return -1;

}

Double_t s::ProperTime(Int_t entry,Int_t loversigma)
{
    //Double_t c=0.2997992458e-4; // speed of light in cm/fs
  Double_t c=0.2997992458e-1; // speed of light in cm/ps
  Double_t p= sqrt(RN_Pxs*RN_Pxs + RN_Pys*RN_Pys + RN_Pzs*RN_Pzs);
  Double_t tau= ((RN_L-loversigma*RN_S)*RN_mass/p)/c; // reduced proper time
  //  Double_t tau= ((RN_Zs-RN_Zp)*RN_mass/p)/c;
   
  return tau;
}

Int_t s::TgtVol(Int_t entry)
{ 
// anal v 2.00  pass2/tgt_pos.ocs and 
// anal v 1.03  pass1/tgt_sz.ocs

  //  Float_t tgt_center_z[5]={6.0889,-4.6519,-3.1552,-1.6611,-0.0538};
  //Float_t tgt_tickness[5]={0.08,0.0572,0.11,0.11,0.11};

  Float_t delta_Z=0.25; // cm
  

 if ((RN_Zs >= -6.1289 - delta_Z && RN_Zs <= (-6.1289 + 0.0800 + delta_Z)) ||
     (RN_Zs >= -4.6233 - delta_Z && RN_Zs <= (-4.6233 + 0.0572 + delta_Z)) ||
     (RN_Zs >= -3.1002 - delta_Z && RN_Zs <= (-3.1002 + 0.11 + delta_Z )) ||
     (RN_Zs >= -1.6061 - delta_Z && RN_Zs <= (-1.6061 + 0.11 + delta_Z )) ||
     (RN_Zs >= -0.1088 - delta_Z && RN_Zs <= (-0.1088 + 0.11 + delta_Z )) )
   {
     return -1;//      cout<<"INSIDE target "<<endl; and we DON'T WANT THAT
   }
 else { 
   return 1; // cout<<"OUTSIDE target "<<endl;
 }
}

Double_t s::Trk_p(Int_t entry, Int_t trk) {
  if (trk==0) {
    Double_t track1_p= sqrt( RN_Pxt1*RN_Pxt1 + RN_Pyt1*RN_Pyt1 + RN_Pzt1*RN_Pzt1 );
    return track1_p;
  }
  if (trk==1) {
    Double_t track2_p= sqrt( RN_Pxt2*RN_Pxt2 + RN_Pyt2*RN_Pyt2 + RN_Pzt2*RN_Pzt2 );
    return track2_p;
  }
  if (trk==2) {
   Double_t  track3_p= sqrt( RN_Pxt3*RN_Pxt3 + RN_Pyt3*RN_Pyt3 + RN_Pzt3*RN_Pzt3 );
    return track3_p;
  }
  if (trk==3) {
   Double_t  track4_p= sqrt( RN_Pxt4*RN_Pxt4 + RN_Pyt4*RN_Pyt4 + RN_Pzt4*RN_Pzt4 );
    return track4_p;
  }
  if (trk==4) {
   Double_t  track5_p= sqrt( RN_Pxt5*RN_Pxt5 + RN_Pyt5*RN_Pyt5 + RN_Pzt5*RN_Pzt5 );
    return track5_p;
  }
  if (trk==5) {
   Double_t  track6_p= sqrt( RN_Pxt6*RN_Pxt6 + RN_Pyt6*RN_Pyt6 + RN_Pzt6*RN_Pzt6 );
    return track6_p;
  }
  if (trk==6) {
  Double_t   track7_p= sqrt( RN_Pxt7*RN_Pxt7 + RN_Pyt7*RN_Pyt7 + RN_Pzt7*RN_Pzt7 );
    return track7_p;
  }
  if (trk==7) {
  Double_t   track8_p= sqrt( RN_Pxt8*RN_Pxt8 + RN_Pyt8*RN_Pyt8 + RN_Pzt8*RN_Pzt8 );
    return track8_p;
  }
  if (trk==8) {
  Double_t   track9_p= sqrt( RN_Pxt9*RN_Pxt9 + RN_Pyt9*RN_Pyt9 + RN_Pzt9*RN_Pzt9 );
    return track9_p;
  }
}

Int_t s::MinOfTrk2Trk3(Int_t entry, float mom_thres) {

Float_t temp1 = s::Trk_p(entry,1);
Float_t temp2 = s::Trk_p(entry,2);
//printf("%f\n",temp1);

if ((temp1>temp2)&&temp2>mom_thres) return 1;
if ((temp2>temp1)&&temp1>mom_thres) return 1;


else return -1;
       
}

Int_t s::MinOfTrk2(Int_t entry, float mom_thres) {
Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, RN_pidt8, 
		  RN_pidt9} ;
Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, RN_mvert_t6, 
		    RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ; 


 Float_t temp1 = s::Trk_p(entry,2);
 // if (temp1>mom_thres && RN_pidt2 != 15 && RN_pidt2 != 24) {OITrk2ID->Fill(RN_pidt2,RN_mvert_t2); return 1; }
 if (temp1>mom_thres) {OITrk2ID->Fill(RN_pidt2,RN_mvert_t2); return 1; }

 else return -1;
 
}

Int_t s::Cut(Int_t entry)
{
// This function may be called from Loop.
// returns  1 if entry is accepted.
// returns -1 otherwise.
   return 1;
}

void s::kp(Int_t entry) {
 Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		   RN_pidt8, RN_pidt9} ;
 Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, RN_mvert_t6, 
		     RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
 Int_t no_of_trk = RN_nt;// e mu pi k p sig Xi ome 
 Float_t temp;          //  8 7  6  5 4 3   2  1 

 for (Int_t i=0; i<no_of_trk; i++) {
   
   if(RN_ID==828) {    

     if (mvert[i]==1 && pid[i]==8 ) { 
       temp=s::Trk_p(entry,i); 
       if(RN_mass>=2.68 && RN_mass<=2.72) {  pi_plus_828->Fill(temp);  }   
       else sb_pi_plus_828->Fill(temp); }
     
     if (mvert[i]==1 && pid[i]==12) { 
       temp=s::Trk_p(entry,i); 
       if(RN_mass>=2.68 && RN_mass<=2.72) { k_minus_828->Fill(temp);  }   
       else sb_k_minus_828->Fill(temp); }
      
     if (mvert[i]==1 && pid[i]==23) {
       temp=s::Trk_p(entry,i); 
       if(RN_mass>=2.68 && RN_mass<=2.72) { X_minus_828->Fill(temp);   } 
       else sb_X_minus_828->Fill(temp); }
   }
   if(RN_ID==829) {     
     
     if (mvert[i]==1 && pid[i]==8) { 
       temp=s::Trk_p(entry,i); 
       if(RN_mass>=2.68 && RN_mass<=2.72) { pi_plus_829->Fill(temp);  } 
       else sb_pi_plus_829->Fill(temp); }
     
     if (mvert[i]==1 && pid[i]==9) {
       temp=s::Trk_p(entry,i); 
       if(RN_mass>=2.68 && RN_mass<=2.72) { pi_minus_829->Fill(temp);   } 
       else sb_pi_minus_829->Fill(temp); }
     
     if (mvert[i]==1 && pid[i]==24) { 
       temp=s::Trk_p(entry,i); 
       if(RN_mass>=2.68 && RN_mass<=2.72) {  Ome_minus_829->Fill(temp);   } 
       else sb_Ome_minus_829->Fill(temp); }
   }
   if(RN_ID==850) {  
     if (mvert[i]==1 && pid[i]==8) { 
       temp=s::Trk_p(entry,i); 
       if(RN_mass>=2.68 && RN_mass<=2.72) {  pi_plus_850->Fill(temp);   } 
       else sb_pi_plus_850->Fill(temp); }

     if (mvert[i]==1 && pid[i]==24) {
       temp=s::Trk_p(entry,i); 
       if(RN_mass>=2.68 && RN_mass<=2.72) {   Ome_minus_850->Fill(temp); }  
       else sb_Ome_minus_850->Fill(temp); }
   }
 }
}

void s::Track_Print(Int_t entry) {
  Float_t m[9] = {RN_masst1, RN_masst2, RN_masst3, RN_masst4, RN_masst5, RN_masst6, RN_masst7, RN_masst8,
		  RN_masst9} ;
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, RN_pidt8,
		    RN_pidt9} ;
  Float_t pid1[9] = {RN_pid1t1, RN_pid1t2, RN_pid1t3, RN_pid1t4, RN_pid1t5, RN_pid1t6, RN_pid1t7, 
		     RN_pid1t8, RN_pid1t9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, RN_mvert_t6, 
		      RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
  Int_t no_of_trk = RN_nt;
  // particle IDs pi+:8 pi-:9 k-:12 p:14 Xi:23
  for(Int_t i=0;i<no_of_trk;i++){
    if (RN_ID==829 && pid[i]==24 && mvert[i]==1) homega->Fill(m[i]);
    //if (RN_ID==828 && pid[i]==23 && mvert[i]==1) hcascade->Fill(m[i]);
  }
}

Int_t s::Angle_Vector(Int_t entry) {
  // e mu pi k prot sig Xi ome  
  // 8 7  6  5 4    3   2  1
  // particle IDs pi+:8 pi-:9 k-:12 p:14 Xi:23
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Float_t pid1[9] = {RN_pid1t1, RN_pid1t2, RN_pid1t3, RN_pid1t4, RN_pid1t5, RN_pid1t6, RN_pid1t7,
		     RN_pid1t8, RN_pid1t9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, RN_mvert_t6,
		      RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ; 
  Int_t no_of_trk = RN_nt;
  Double_t vec_x[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9};
  Double_t vec_y[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9};
  Double_t vec_z[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9};

  Int_t kaon_index=-1; Int_t pion_index=-1; Int_t Xi_index=-1;
  
  for(Int_t i=0; i<no_of_trk ; i++) {
    if (pid[i]==12 && mvert[i]==1) kaon_index=i; // found kaon track ind
    if (pid[i]==23 && mvert[i]==1) Xi_index=i; // found cascade track ind
    //    cout<<"kaon_index "<<kaon_index<<" Xi_index "<<Xi_index<<endl;    
  }

  if(kaon_index>=0 && Xi_index>=0) {

TVector3 v1(vec_x[kaon_index], vec_y[kaon_index], vec_z[kaon_index]); //kaon momentum vector
TVector3 v2(vec_x[Xi_index], vec_y[Xi_index], vec_z[Xi_index]); // Xi momentum vector
 Double_t ang = (180./TMath::Pi()) * v1.Angle(v2);
 if(ang<=0.8) return 1;
 else return -1;
  }
}

Int_t s::Event_Reject(Int_t entry) {
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Float_t pid1[9] = {RN_pid1t1, RN_pid1t2, RN_pid1t3, RN_pid1t4, RN_pid1t5, RN_pid1t6, RN_pid1t7, 
		     RN_pid1t8, RN_pid1t9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, RN_mvert_t6, 
		      RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
  Int_t no_of_trk = RN_nt;
  // e mu pi k prot sig Xi ome  
  // 8 7  6  5 4    3   2  1  // particle IDs pi+:8 pi-:9 k-:12 p:14
  for(Int_t i=0;i<RN_nt;i++){
    //    cout<<RN_run<<" "<<RN_event<<" "<<RN_mass<<" "<<RN_ID<<endl;   
  } 
}


Int_t s::Mvert_Part(Int_t entry, Int_t ver_indx, Float_t Momentum, const char *part_name) {
  Int_t part_index;  Int_t no_of_trk = RN_nt; Int_t part_id;
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Float_t pid1[9] = {RN_pid1t1, RN_pid1t2, RN_pid1t3, RN_pid1t4, RN_pid1t5, RN_pid1t6, 
		     RN_pid1t7, RN_pid1t8, RN_pid1t9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
 Int_t part_indx[3] = {-1,-1,-1};  Int_t j=0;
 if (part_name=="k-")  part_id=12;      // k- =12, 
 if (part_name=="pi+") part_id=8;      //pi+ =8,
 if (part_name=="pi-") part_id=9;     // pi- =9, 
 if (part_name=="p")   part_id=14;   //p  =14
 if (part_name=="ome-") part_id=24; // ome-  =24
  
 for(Int_t i=0; i<no_of_trk ; i++) {
   if (pid[i]==part_id && mvert[i]==ver_indx) {part_indx[j]=i; j++;}
  
 } 
 if(s::Trk_p(entry,part_indx[0])>Momentum) return 1;
 else
   return -1;
 
 /* this part is "old",  in fact bugged. i modified it as in above!!! 11/18/02 
    for(Int_t i=0; i<no_of_trk; i++) {
    if (pid[i]==part_id) {
    //     cout<<part_name<<" "<<part_id<<" "<<mvert[i]<<endl;
    part_index=i; // found requested particle's track index
    break;}
    }
    //now check if this particle is coming from requested mvert ?? 
    if (mvert[part_index]==ver_indx) {   // yes , it does...so ...
    if( s::Trk_p(entry,part_index+1)>=Momentum)  return 1; //momentum threshold
    else 
    return -1;
    }
 */
}

void s::Lambda0(Int_t entry) { //Calculate the Lambda0 mass from proton and pi- at mvert=2
   Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
                     RN_pidt8, RN_pidt9} ;
   Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
                       RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
   Float_t mass[9] = {RN_masst1, RN_masst2, RN_masst3, RN_masst4, RN_masst5, RN_masst6, 
                      RN_masst7, RN_masst8, RN_masst9} ;
   Int_t no_of_trk = RN_nt;
   Double_t p_x[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9};
   Double_t p_y[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9};
   Double_t p_z[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9};
  
  Double_t mass_proton=0.93827231; 
  Double_t mass_pi_minus=0.1395679; 
  Double_t E_proton; Double_t E_pi_minus; //Energies
  Double_t P_proton; Double_t P_pi_minus; //Momenta
  Int_t pro_indx=-1; Int_t pi_minus_indx=-1; //Both are decay products of Lambda_0 from mvert=2
  
  for(Int_t i=0; i<no_of_trk; i++) {// particle IDs pi+:8 pi-:9 k-:12 p:14
    if(pid[i]==9 && mvert[i]==3) pi_minus_indx=i; //found the pi_minus from Lambda0 decay
    if(pid[i]==14 ) pro_indx=i; //found the proton from Lambda0 decay
  }

  if(pro_indx != -1 && pi_minus_indx != -1) { //if both exist => calculate kinematics variables

P_proton=sqrt(pow(p_x[pro_indx],2)+ pow(p_y[pro_indx],2)+ pow(p_z[pro_indx],2));
E_proton=sqrt(pow(mass[pro_indx],2) + pow(P_proton,2));
P_pi_minus=sqrt(pow(p_x[pi_minus_indx],2)+ pow(p_y[pi_minus_indx],2)+ pow(p_z[pi_minus_indx],2));
E_pi_minus=sqrt(pow(mass[pi_minus_indx],2) + pow(P_pi_minus,2));


 Double_t m1=mass_proton;  // 1 : proton 
 Double_t m2=mass_pi_minus;// 2 : pion
 Double_t E1=E_proton;
 Double_t E2=E_pi_minus;
 Double_t p1x=p_x[pro_indx]; Double_t p1y=p_y[pro_indx]; Double_t p1z=p_z[pro_indx];
 Double_t p2x=p_x[pi_minus_indx]; Double_t p2y=p_y[pi_minus_indx]; Double_t p2z=p_z[pi_minus_indx
];

 Double_t M = sqrt (pow(m1+m2,2) + 2*(E1*E2 - m1*m2 - (p1x*p2x + p1y*p2y + p1z*p2z) ) );
 
hlambda->Fill(M);
  }

}

  // particle IDs pi+:8 pi-:9 k-:12 p:14 Ome- =24 X-=23  sigma + ==19 , sigma - = 21
Double_t s::Cos_Theta(Int_t entry) {
  Float_t m[9] = {RN_masst1,RN_masst2,RN_masst3,RN_masst4,RN_masst5,RN_masst6,RN_masst7,
		  RN_masst8,RN_masst9} ;
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		     RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
  Float_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
  Float_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ;
  Float_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9} ;
  Float_t cos_theta; Int_t no_of_trk = RN_nt; Int_t i,j,k;
  Float_t e1; Float_t e1_cms;  Float_t e2; Float_t e2_cms; Float_t mA;
  Float_t p1; Float_t p1_cms;  Float_t pA; Float_t pA_cms; Float_t eA;
  Int_t pi[3] ={-1,-1,-1}; Int_t ome[1] ={-1};
  Float_t m1,m2;
  for(Int_t i=0; i<no_of_trk ; i++) {
    if (pid[i]==8 && mvert[i]==1) {pi[j]=i; j++;}
    else 
      if (pid[i]==24 && mvert[i]==1) {ome[k]=i;}
  }   
mA= RN_mass; 
m1= m[pi[0]];
m2= m[ome[0]];
p1= sqrt(pow(px[pi[0]],2) + pow(py[pi[0]],2) + pow(pz[pi[0]],2));
p1_cms= sqrt(pow(mA,2) - pow((m1+m2),2))  * sqrt(pow(mA,2) - pow((m1-m2),2)) / (2 * mA);
e1_cms= (pow(mA,2) + pow(m1,2) - pow(m2,2))/(2*mA);
e1 = sqrt(m1*m1 + p1*p1);
pA = sqrt(pow(RN_Pxs,2) + pow(RN_Pys,2) + pow(RN_Pzs,2));
eA = sqrt(mA*mA + pA*pA);

cos_theta= (e1*mA - e1_cms*eA)/(p1_cms*pA);

 if((RN_mass>=2.685 && RN_mass<=2.715)) hcos_teti->Fill(cos_theta);
 if(!(RN_mass>=2.685 && RN_mass<=2.715)) hcos_teto->Fill(cos_theta);
 return cos_theta;
}



void s::Slice2DHisto(TH2F *hp) {
   gStyle->SetLabelSize(0.06,"x");
   gStyle->SetLabelSize(0.06,"y");
   gStyle->SetFrameFillColor(38);
   gStyle->SetTitleW(0.6);
   gStyle->SetTitleH(0.1);
// Connect the input file and get the 2-d histogram in memory
   TH2F *hpxpy =(TH2F*)hp->Clone();
   hpxpy->SetName("hpxpy");
   // Create a canvas and divide it
   TCanvas *c7 = new TCanvas("c7","c7",700,500);
   c7->SetFillColor(42);
   c7->Divide(2,1);
   c7->cd(1);
   TPad *left = (TPad*)gPad;
   left->Divide(1,2);
// Draw 2-d original histogram
   left->cd(1);
   gPad->SetTopMargin(0.12);
   gPad->SetFillColor(33);
   hpxpy->Draw("box");
   hpxpy->GetXaxis()->SetLabelSize(0.04);
   hpxpy->GetYaxis()->SetLabelSize(0.04);
   hpxpy->SetMarkerColor(kYellow);
// Fit slices projected along Y fron bins in X [7,32]
TF1 *pullfcn = new TF1("pullfcn","gaus",-10,10);    // Create fitting function 
pullfcn->SetRange(-8.,8.);
hpxpy->FitSlicesY(pullfcn,0,15,0);
// Show fitted "mean" for each slice
   left->cd(2);
   gPad->SetFillColor(33);
   hpxpy_0->Draw("box");
   c7->cd(2);
   TPad *right = (TPad*)gPad;
   right->Divide(1,2);
   right->cd(1);
   gPad->SetTopMargin(0.12);
   gPad->SetLeftMargin(0.15);
   gPad->SetFillColor(33);
   hpxpy_1->Draw();
// Show fitted "sigma" for each slice
   right->cd(2);
   gPad->SetTopMargin(0.12);
   gPad->SetLeftMargin(0.15);
   gPad->SetFillColor(33);
   // hpxpy_2->SetMinimum(0.);
   hpxpy_2->Draw();
//attributes
   hpxpy_0->SetLineColor(kYellow);
   hpxpy_1->SetLineColor(kYellow);
   hpxpy_2->SetLineColor(kYellow);
   hpxpy_0->SetMarkerColor(kRed);
   hpxpy_1->SetMarkerColor(kRed);
   hpxpy_2->SetMarkerColor(kRed);
   hpxpy_0->SetMarkerStyle(21);
   hpxpy_1->SetMarkerStyle(21);
   hpxpy_2->SetMarkerStyle(21);
   hpxpy_0->SetMarkerSize(0.6);
   hpxpy_1->SetMarkerSize(0.6);
   hpxpy_2->SetMarkerSize(0.6);
}

void s::Histo_Normalizer() {
  Double_t norm;
  
  norm= hAngle_in->GetEntries(); if(norm) hAngle_in->Scale(1/norm);
  norm= hAngle_out->GetEntries(); if(norm) hAngle_out->Scale(1/norm);
  
  norm= hSumpt2_in->GetEntries(); if(norm) hSumpt2_in->Scale(1/norm);
  norm= hSumpt2_out->GetEntries(); if(norm) hSumpt2_out->Scale(1/norm);
  
  norm= hOIpt_in->GetEntries(); if(norm) hOIpt_in->Scale(1/norm);
  norm= hOIpt_out->GetEntries(); if(norm) hOIpt_out->Scale(1/norm); 
}

void s::ProperTimePlot() {
TCanvas *c2 = new TCanvas("c2","Proper Time distrutions",0,0,600,600);   //canvas 2
 c2->Divide(2,2);
 c2->cd(1);h827ProperTime->Draw("e");
 c2->cd(2);h828ProperTime->Draw("e");
 c2->cd(3);h829ProperTime->Draw("e");
 c2->cd(4);h850ProperTime->Draw("e");


}

void s::OmeZs2ProTime() {
TCanvas *c7 = new TCanvas("c7","#Omega^{-} ",0,0,600,600);   //canvas 6
TPostScript *ps4 = new TPostScript("MC_time.ps",111);
c7->Clear();
c7->Divide(2,2);
//c7->Update();ps4->NewPage();
c7->cd(1);c7_1->SetLogy(1);ht829->Draw();
c7->cd(2);c7_2->SetLogy(1);ht1829->Draw();
c7->cd(3);c7_3->SetLogy(1);ht2829->Draw(); 
c7->Update();ps4->NewPage();

//c7->cd(1);c7_1->SetLogy(1);ht3829->Draw();
//c7->cd(2);c7_2->SetLogy(1);ht3850->Draw();
//c7->cd(3);c7_3->SetLogy(1);ht4829->Draw(); 
//c7->cd(4);c7_4->SetLogy(1);ht4850->Draw(); 
c7->cd(1);c7_1->SetLogy(1);ht22829->Draw();
c7->cd(2);c7_2->SetLogy(1);ht22850->Draw();
c7->cd(3);c7_3->SetLogy(1);ht222829->Draw(); 
c7->cd(4);c7_4->SetLogy(1);ht222850->Draw(); 
c7->Update();ps4->NewPage();


 ht829->Clone("ht829norm"); ht1829->Clone("ht1829norm"); ht2829->Clone("ht2829norm");
 norm=ht829norm->GetEntries();ht829norm->Scale(1/norm);
 norm=ht1829norm->GetEntries();ht1829norm->Scale(1/norm);
 norm=ht2829norm->GetEntries();ht2829norm->Scale(1/norm);

c7->cd(1);c7_1->SetLogy(1);ht829norm->Draw();
c7->cd(2);c7_2->SetLogy(1);ht1829norm->Draw();
c7->cd(3);c7_3->SetLogy(1);ht2829norm->Draw(); c7->Update();ps4->NewPage();

 TF1 *fun1 = new TF1("fun1","exp(-x/.82)",0,5.); // Oiii MC
 TH1F *hfun1 =new TH1F("hfun1", "Ome- production propertime",25,0,5.0);
 TH1F *hfun2 =new TH1F("hfun2", "Ome- production propertime",25,0,5.0);
 TH1F *correctionfun1 =new TH1F("correctionfun1", "exp/Ome- all",25,0,5.);
 TH1F *correctionfun2 =new TH1F("correctionfun2", "exp/ome- 100<z<200",25,0,5.);
 correctionfun1->Sumw2();correctionfun2->Sumw2();
 //hfun1->Eval(fun1);hfun1->Scale(1/ht1829->GetEntries());
 //hfun2->Eval(fun1);hfun2->Scale(1/ht2829->GetEntries());
 hfun1->FillRandom("fun1",1000000); //
 hfun2->FillRandom("fun1",1000000); // create a perfect exponential 
 hfun1->Scale(ht1829->GetEntries()/1000000);// scale with number of entry
 hfun2->Scale(ht2829->GetEntries()/1000000);

 correctionfun1->Divide(hfun1,ht1829,1,1); 
 correctionfun2->Divide(hfun2,ht2829,1,1); 
 c7->cd(1);correctionfun1->Draw("e");
 c7->cd(2);correctionfun2->Draw("e");

c7->Update();ps4->NewPage();

 hfun1->FillRandom("fun1",1000000); //
 hfun2->FillRandom("fun1",1000000); // create a perfect exponential 
 TH1F *correctionfun3 =new TH1F("correctionfun3", "829 exp/Ome- signalband",25,0,5.);
 TH1F *correctionfun4 =new TH1F("correctionfun4", "850 exp/ome- signalband ",25,0,5.);
 hfun1->Scale(ht829sbs->GetEntries()/1000000);// scale with number of entry
 hfun2->Scale(ht850sbs->GetEntries()/1000000);

 norm1=ht829gen->GetEntries();  norm2=ht22829->GetEntries(); 
 ht829gen->Scale(norm2/norm1);
 
 norm1=ht850gen->GetEntries();  norm2=ht22850->GetEntries(); 
 ht850gen->Scale(norm2/norm1);

 correctionfun3->Divide(ht829gen,ht22829,1,1); 
 correctionfun4->Divide(ht850gen,ht22850,1,1); 

 c7->cd(1);correctionfun3->Draw("e");
 c7->cd(2); correctionfun4->Draw("e");

 ps4->Close();
}



Int_t s::pid829(Int_t entry) {
  Int_t no_of_trk = RN_nt;
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Int_t pid1[9] = {RN_pid1t1, RN_pid1t2, RN_pid1t3, RN_pid1t4, RN_pid1t5, RN_pid1t6, RN_pid1t7,
		   RN_pid1t8, RN_pid1t9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
  Double_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
  Double_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ; 
  Double_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9} ;
  Int_t kaon_minus= -1;  Int_t pion_minus= -1;Int_t pion_plus[2]= {-1,-1}; Int_t Ome_minus = -1; 
  Int_t j=0; Int_t pion[3]={-1,-1,-1};
  // particle IDs pi+:8 pi-:9 k-:12 p:14 Ome_minus=24 X_minus=23
  for(Int_t i=0; i<RN_nt ; i++) {
    if( mvert[i]==1 && (pid[i]==8 || pid[i]==9)) 
      {
	pion[j]=i;j++;
      }
  }
  
  Double_t pion1_mom = Trk_p(entry,pion[0]);
  Double_t pion2_mom = Trk_p(entry,pion[1]);
  Double_t pion3_mom = Trk_p(entry,pion[2]);

  //if(pion1_mom>12. && pion2_mom>12. && pion3_mom>12.) return 1;
  if(pion1_mom>8. && pion2_mom>8. && pion3_mom>8.) return 1;
  else return -1;
  /* Int_t n=3; // Bubble sort the momentum in ascending order
     Double_t a[] = {pion1_mom,pion2_mom,pion3_mom};
     for(i=0; i<n; i++){
     Bubble(a,n);
     }
    if(a[0]>8. && a[1]>12. && a[2]>15.) return 1;
    else return -1;
  */
}
Int_t s::pid850(Int_t entry) {
  Int_t no_of_trk = RN_nt;
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Int_t pid1[9] = {RN_pid1t1, RN_pid1t2, RN_pid1t3, RN_pid1t4, RN_pid1t5, RN_pid1t6, RN_pid1t7,
		     RN_pid1t8, RN_pid1t9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
  Double_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
  Double_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ; 
  Double_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9} ;
  Int_t kaon_minus= -1;  Int_t pion_minus= -1; Int_t pion_plus= -1; Int_t Ome_minus = -1; Int_t j=0;
  Int_t pion[3]={-1,-1,-1}; Int_t flag=-1;
  // particle IDs pi+:8 pi-:9 k-:12 p:14 Ome_minus=24 X_minus=23
   for(Int_t i=0; i<no_of_trk ; i++) {
     if(mvert[i]==1 && (pid[i]==8 || pid[i]==9) )  {
       Double_t momentum= Trk_p(entry,i); 
       //Double_t momentum= sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]); 
       
     }
   }
   if(momentum>12.) return 1;
   else return -1;
}

Int_t s::pid828(Int_t entry) {
  Int_t no_of_trk = RN_nt;
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Int_t pid1[9] = {RN_pid1t1, RN_pid1t2, RN_pid1t3, RN_pid1t4, RN_pid1t5, RN_pid1t6, RN_pid1t7,
		     RN_pid1t8, RN_pid1t9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
  Double_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
  Double_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ; 
  Double_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9} ;
  Int_t kaon_minus= -1;  Int_t pion_minus= -1; Int_t pion_plus= -1; Int_t Ome_minus = -1; Int_t j=0;
  Int_t pion[3]={-1,-1,-1}; Int_t flag=-1;
  // particle IDs pi+:8 pi-:9 k-:12 p:14 Ome_minus=24 X_minus=23
   for(Int_t i=0; i<no_of_trk ; i++) {
     if(mvert[i]==1 && pid[i]==8)  {
       pion[j]=i;
       j++;
     }
   }
   Double_t pimom1= Trk_p(entry,pion[0]);
   Double_t pimom2= Trk_p(entry,pion[1]);
   
   if(pimom1>8.&& pimom2>8.) return 1;
   else return -1;
}

Double_t s::Oi_p_asym(Int_t entry){
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ; // particle IDs pi+:8 pi-:9 k-:12 p:14 Ome_minus=24 X_minus=23
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
  Int_t pion=-1;Int_t ome=-1;
  for(Int_t i=0; i<RN_nt ; i++) {
    if(mvert[i]==1 && pid[i]==8) pion=i;
    if(mvert[i]==1 && pid[i]==24) ome=i;
  }
  Double_t ppi=Trk_p(entry,pion);
  Double_t omepi=Trk_p(entry,ome);
  Double_t p_asym=(ppi-omepi)/(ppi+omepi);//momentum asymmetry
  //cout<<p_asym<<endl;
  return  p_asym;
}


Int_t s::kk827(Int_t entry) {
 
  Float_t m[9] = {RN_masst1,RN_masst2,RN_masst3,RN_masst4,RN_masst5,RN_masst6,RN_masst7,
		  RN_masst8,RN_masst9} ;
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Float_t pid1[9] = {RN_pid1t1, RN_pid1t2, RN_pid1t3, RN_pid1t4, RN_pid1t5, RN_pid1t6, RN_pid1t7,
		     RN_pid1t8, RN_pid1t9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
 Int_t no_of_trk = RN_nt; Int_t j=0,k=0,l=0;
 Int_t kaon[3]={-1,-1,-1};
 // particle IDs pi+:8 pi-:9 k-:12 p:14 Ome_minus=24 X_minus=23
  // e mu pi k prot sig Xi ome  
  // 8 7  6  5 4    3   2  1
  for(Int_t i=0; i<no_of_trk ; i++) {
    if (pid[i]==12 && mvert[i]==1) {kaon[j]=i; j++;}
    // int((int(pid1[part_index])%100000) / 10000)   ; //kaon likelihood
  }
  //  cout<<int((int(pid1[kaon[0]])%100000) / 10000)<<" "<<  int((int(pid1[kaon[1]])%100000) / 10000)   <<endl;
 
  if(int((int(pid1[kaon[0]])%100000) / 10000)>5 && int((int(pid1[kaon[1]])%100000) / 10000)>5) return 1;
  else 
    return -1;
  
}


Double_t s::oiiipt(Int_t entry) {//829 
 Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, RN_pidt8, RN_pidt9} ;
 Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		     RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
 Float_t m[9] = {RN_masst1,RN_masst2,RN_masst3,RN_masst4,RN_masst5,RN_masst6,RN_masst7,RN_masst8,RN_masst9} ;
 Float_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
 Float_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ;
 Float_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt9,RN_Pzt9} ;
 Int_t pi[3] = {-1,-1,-1};  Int_t pi_mi[3] = {-1,-1,-1}; Int_t ome_minus = -1; 
 Int_t j=0,k=0;

 for(Int_t i=0; i<RN_nt ; i++) {
   if ((pid[i]==8||pid[i]==9) && mvert[i]==1) {pi[j]=i; j++;}
 }
 //cout<<entry<<" "<<pi[0]<<pi[1]<<pi[2]<<endl;

 double p1x=px[pi[0]]; double p1y=py[pi[0]]; double p1z=pz[pi[0]];
 double p2x=px[pi[1]]; double p2y=py[pi[1]]; double p2z=pz[pi[1]];
 double p3x=px[pi[2]]; double p3y=py[pi[2]]; double p3z=pz[pi[0]];
 
 // a_vector : omegac0 momentum
 // p_vector : Sum of 3 pi momentum
 // bt       : transverse component of b_vector (w.r.t p_vector) 
 Double_t ax=RN_Pxs;  Double_t ay=RN_Pys;  Double_t az=RN_Pzs;
 Double_t a =sqrt(ax*ax + ay*ay + az*az); //length of a
 
 Double_t psumx=p1x+p2x+p3x;Double_t psumy=p1y+p2y+p3y;Double_t psumz=p1z+p2z+p3z; // p = p1+p2+p3
 Double_t psum=sqrt(psumx*psumx + psumy*psumy + psumz*psumz); //length of p

 Double_t adotp = ax*psumx + ay*psumy + az*psumz;
 Double_t bt    = sqrt( psum*psum - pow(adotp/a,2)    );

//Double_t bt   =sqrt(psumx*psumx + psumy*psumy);
 //cout<<bt<<endl;
 return bt;
}

Double_t s::oipt(Int_t entry) { //To Calculate Pt of Pi+   Nov.,17,2002
 // particle IDs pi+:8 pi-:9 Ome_minus=24
 Int_t no_of_trk = RN_nt;
 Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, RN_pidt8, RN_pidt9} ;
 Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		     RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
 Float_t m[9] = {RN_masst1,RN_masst2,RN_masst3,RN_masst4,RN_masst5,RN_masst6,RN_masst7,RN_masst8,RN_masst9} ;
 Float_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
 Float_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ;
 Float_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt9,RN_Pzt9} ;
 Int_t pi_pl[3] = {-1,-1,-1};  Int_t ome_min = -1;  Int_t j=0,k=0;
 //Double_t bt, adotb, eta, ax,ay,az, bx,by,bz;
 for(Int_t i=0; i<no_of_trk ; i++) {
   if (pid[i]==8 && mvert[i]==1) {pi_pl[j]=i; j++;}
     if (pid[i]==24 && mvert[i]==1) {ome_min=i; }
 } 
 //cout<<pi_pl[0]<<" "<<pi_pl[1]<<" "<<pi_pl[2]<<" "<<endl;

 // a_vector : omegac0 momentum
 // b_vector : pi+     momentum
 // bt       : transverse component of b_vector (w.r.t a_vector) 
 Double_t ax=RN_Pxs;  Double_t ay=RN_Pys;  Double_t az=RN_Pzs;
 Double_t a =sqrt(ax*ax + ay*ay + az*az); //length of a
 
 Double_t bx=px[pi_pl[0]]; Double_t by=py[pi_pl[0]]; Double_t bz=pz[pi_pl[0]];
 Double_t b=sqrt(bx*bx + by*by + bz*bz); //length of b

 Double_t adotb = ax*bx + ay*by + az*bz;

 // Double_t  eta   = adotb / (ax*ax + ay*ay + az*az);
 // Double_t bt    = sqrt((bx-eta*ax)*(bx-eta*ax) + (by-eta*ay)*(by-eta*ay) +(bz-eta*az)*(bz-eta*az));

 Double_t bt    = sqrt( b*b - pow(adotb/a,2) );
 // cout<<bt<<endl;
 return bt;

}

Double_t s::sum_pt2(Int_t entry) {
  //sum of pt**2 wrt Omegac0 momentum
 Int_t no_of_trk = RN_nt;
 Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, RN_pidt8, RN_pidt9} ;
 Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		     RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
 Float_t m[9] = {RN_masst1,RN_masst2,RN_masst3,RN_masst4,RN_masst5,RN_masst6,RN_masst7,RN_masst8,RN_masst9} ;
 Float_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
 Float_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ;
 Float_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt9,RN_Pzt9} ;
 Int_t pi[3] = {-1,-1,-1};  Int_t ome = -1; 
 Int_t j=0,k=0;
  for(Int_t i=0; i<no_of_trk ; i++) {
   if ((pid[i]==8||pid[i]==9) && mvert[i]==1) {pi[j]=i; j++;}
     if (pid[i]==24 && mvert[i]==1) {ome=i; }
 }
 Double_t p  = sqrt(RN_Pxs*RN_Pxs       + RN_Pys*RN_Pys       + RN_Pzs*RN_Pzs       );
 Double_t p1 = sqrt(px[pi[0]]*px[pi[0]] + py[pi[0]]*py[pi[0]] + pz[pi[0]]*pz[pi[0]] );
 Double_t p2 = sqrt(px[pi[1]]*px[pi[1]] + py[pi[1]]*py[pi[1]] + pz[pi[1]]*pz[pi[1]] );
 Double_t p3 = sqrt(px[pi[2]]*px[pi[2]] + py[pi[2]]*py[pi[2]] + pz[pi[2]]*pz[pi[2]] );
 Double_t p4 = sqrt(px[ome]*px[ome]     + py[ome]*py[ome]     + pz[ome]*pz[ome]     );

 Double_t p_dot_p1 = px[pi[0]]*RN_Pxs + py[pi[0]]*RN_Pys + pz[pi[0]]*RN_Pzs;
 Double_t p_dot_p2 = px[pi[1]]*RN_Pxs + py[pi[1]]*RN_Pys + pz[pi[1]]*RN_Pzs;
 Double_t p_dot_p3 = px[pi[2]]*RN_Pxs + py[pi[2]]*RN_Pys + pz[pi[2]]*RN_Pzs;
 Double_t p_dot_p4 = px[ome]  *RN_Pxs + py[ome]  *RN_Pys + pz[ome]  *RN_Pzs;

 Double_t pt1= sqrt(p1*p1 - pow(p_dot_p1/p,2));
 Double_t pt2= sqrt(p2*p2 - pow(p_dot_p2/p,2));
 Double_t pt3= sqrt(p3*p3 - pow(p_dot_p3/p,2));
 Double_t pt4= sqrt(p4*p4 - pow(p_dot_p4/p,2));

 //cout<<p1<<" "<<p2<<" "<<p3<<" "<<p4<<" "<<endl;
 //cout<<pt1<<" "<<pt2<<" "<<pt3<<" "<<pt4<<" "<<endl;

 // Double_t s_pt2= p1*p1 + p2*p2 + p3*p3 + p4*p4 -                      //sum of pt2
 // (  pow(px[pi[0]]*RN_Pxs +py[pi[0]]*RN_Pys +pz[pi[0]]*RN_Pzs,2)+
 //    pow(px[pi[1]]*RN_Pxs +py[pi[1]]*RN_Pys +pz[pi[1]]*RN_Pzs,2)+
 //    pow(px[pi[2]]*RN_Pxs +py[pi[2]]*RN_Pys +pz[pi[2]]*RN_Pzs,2)+
 //    pow(px[ome]*RN_Pxs +py[ome]*RN_Pys +pz[ome]*RN_Pzs,2))/pow(p,2);

  Double_t s_pt2= pt1*pt1 + pt2*pt2 + pt3*pt3 + pt4*pt4 ;            //sum of pt2 
  //cout<<"s_pt2 "<<s_pt2<<endl;
 return s_pt2;

}
Double_t s::sum_pt(Int_t entry) {
  //sum of pt of 3 pions wrt Omegac0 momentum
 Int_t no_of_trk = RN_nt;
 Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, RN_pidt8, RN_pidt9} ;
 Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		     RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
 Float_t m[9] = {RN_masst1,RN_masst2,RN_masst3,RN_masst4,RN_masst5,RN_masst6,RN_masst7,RN_masst8,RN_masst9} ;
 Float_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
 Float_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ;
 Float_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt9,RN_Pzt9} ;
Int_t pi[3] = {-1,-1,-1};  Int_t ome = -1; 
 Int_t j=0,k=0;
  for(Int_t i=0; i<no_of_trk ; i++) {
   if ((pid[i]==8||pid[i]==9) && mvert[i]==1) {pi[j]=i; j++;}
     if (pid[i]==24 && mvert[i]==1) {ome=i; }
 }
 Double_t p  = sqrt(RN_Pxs*RN_Pxs       + RN_Pys*RN_Pys       + RN_Pzs*RN_Pzs) ;
 Double_t p1 = sqrt(px[pi[0]]*px[pi[0]] + py[pi[0]]*py[pi[0]] + pz[pi[0]]*pz[pi[0]] );
 Double_t p2 = sqrt(px[pi[1]]*px[pi[1]] + py[pi[1]]*py[pi[1]] + pz[pi[1]]*pz[pi[1]] );
 Double_t p3 = sqrt(px[pi[2]]*px[pi[2]] + py[pi[2]]*py[pi[2]] + pz[pi[2]]*pz[pi[2]] );
 Double_t p4 = sqrt(px[ome]*px[ome]     +  py[ome]*py[ome]    + pz[ome]*pz[ome] );

 Double_t p_dot_p1 = px[pi[0]]*RN_Pxs +py[pi[0]]*RN_Pys +pz[pi[0]]*RN_Pzs;
 Double_t p_dot_p2 = px[pi[1]]*RN_Pxs +py[pi[1]]*RN_Pys +pz[pi[1]]*RN_Pzs;
 Double_t p_dot_p3 = px[pi[2]]*RN_Pxs +py[pi[2]]*RN_Pys +pz[pi[2]]*RN_Pzs;
 Double_t p_dot_p4 = px[ome]*RN_Pxs   +py[ome]*RN_Pys   +pz[ome]*RN_Pzs;

 Double_t pt1= sqrt(p1*p1 - pow(p_dot_p1/p,2));
 Double_t pt2= sqrt(p2*p2 - pow(p_dot_p2/p,2));
 Double_t pt3= sqrt(p3*p3 - pow(p_dot_p3/p,2));
 Double_t pt4= sqrt(p4*p4 - pow(p_dot_p4/p,2));

 // Double_t spt=pt1+pt2+pt3;
 Double_t spt=sqrt((px[pi[0]]+px[pi[1]]+px[pi[2]])**2 +(py[pi[0]]+py[pi[1]]+py[pi[2]])**2 ); //sum of pt of 3 pions 
return spt;

}

Double_t s::pts(Int_t entry) {
  // pt along the line of flight
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, RN_pidt8, RN_pidt9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
 Float_t m[9] = {RN_masst1,RN_masst2,RN_masst3,RN_masst4,RN_masst5,RN_masst6,RN_masst7,RN_masst8,RN_masst9} ;
 Float_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
 Float_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ;
 Float_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt9,RN_Pzt9} ;
 Int_t pi[3] = {-1,-1,-1};  Int_t ome = -1; 

 Double_t lx = RN_Xs2 - RN_Xs;
 Double_t ly = RN_Ys2 - RN_Ys;
 Double_t lz = RN_Zs2 - RN_Zs;
 Double_t ls = sqrt(lx*lx + ly*ly + lz*lz);
 Double_t pdotls=RN_Pxs*lx + RN_Pys*ly + RN_Pzs*lz ;

 Double_t ptx=RN_Pxs-lx*pdotls/pow(ls,2);
 Double_t pty=RN_Pys-ly*pdotls/pow(ls,2);

  return sqrt( ptx*ptx + pty*pty );

}




Double_t s::ReCalc4(Int_t entry) { // written for 828

  Double_t p_x[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9};
  Double_t p_y[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9};
  Double_t p_z[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9};
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
  Double_t m[9] = {RN_masst1,RN_masst2,RN_masst3,RN_masst4,RN_masst5,RN_masst6,RN_masst7,RN_masst8,RN_masst9};
  Double_t kaonmass=0.493677; Double_t pionmass=0.1395679; Double_t sigmamass=1.197449; Double_t Ximass=1.32131;
  Double_t p1,p2,p3,p4,E1,E2,E3,E4;
  Int_t piplus[2]={-1,-1}; Int_t kminus=-1; Int_t Xi=-1;Int_t piminus=-1; Int_t Ome=-1;
  Int_t j=0,k=0;
  for(Int_t i=0; i<RN_nt; i++) {// particle IDs pi+:8 pi-:9 k-:12 p:14
    if(mvert[i]==1) {
      if(pid[i]==8) {piplus[j]=i; j++;}
      if(pid[i]==9) piminus=i; //cout<<"k"<<kminus<<endl;
      if(pid[i]==12) kminus=i; //cout<<"k"<<kminus<<endl;
      if(pid[i]==23) Xi=i; 
      if(pid[i]==24) Ome=i; 
    }
  } //must have found pion+'s and kaon from vertex 2
  //now ReCalculate the Omegac0 mass with K- replaced with Pi-
  p1=sqrt(pow(p_x[piplus[0]],2) + pow(p_y[piplus[0]],2) + pow(p_z[piplus[0]],2)); // 1st pi+
  p2=sqrt(pow(p_x[piplus[1]],2) + pow(p_y[piplus[1]],2) + pow(p_z[piplus[1]],2)); // 2nd pi+
  //  p3=sqrt(pow(p_x[kminus],2) + pow(p_y[kminus],2) + pow(p_z[kminus],2));          // K-
  p3=sqrt(pow(p_x[piminus],2) + pow(p_y[piminus],2) + pow(p_z[piminus],2));          // pi-
  //  p4=sqrt(pow(p_x[Xi],2) + pow(p_y[Xi],2) + pow(p_z[Xi],2)); // Xi- from vertex information
  p4=sqrt(pow(p_x[Ome],2) + pow(p_y[Ome],2) + pow(p_z[Ome],2)); // Ome- from vertex information
  
  E1=sqrt(pow(m[piplus[0]],2) + pow(p1,2));
  E2=sqrt(pow(m[piplus[1]],2) + pow(p2,2));
  E3=sqrt(pow(pionmass,2) + pow(p3,2)); // 
  //  E3=sqrt(pow(m[kminus],2) + pow(p3,2)); 
  E4=sqrt(pow(Ximass,2) + pow(p4,2)); //Xi energy
  // E4=sqrt(pow(sigmamass,2) + pow(p4,2)); // Replace Xi- with Ome-
  
  Double_t newmass= sqrt(pow(E1+E2+E3+E4,2) - 
			 pow(p_x[piplus[0]]+p_x[piplus[1]]+p_x[piminus]+p_x[Ome],2) -
			 pow(p_y[piplus[0]]+p_y[piplus[1]]+p_y[piminus]+p_y[Ome],2) -
			 pow(p_z[piplus[0]]+p_z[piplus[1]]+p_z[piminus]+p_z[Ome],2)
			 );
    return newmass;
  //cout<<"newmass"<<newmass<<endl;
}
Double_t s::ReCalc2(Int_t entry) { // written for 850

  Double_t p_x[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9};
  Double_t p_y[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9};
  Double_t p_z[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9};
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
  Double_t m[9] = {RN_masst1,RN_masst2,RN_masst3,RN_masst4,RN_masst5,RN_masst6,RN_masst7,RN_masst8,RN_masst9};
  Double_t kaonmass=0.493677; Double_t pionmass=0.13957018; 
  Double_t Ximass=1.32131;    Double_t sigmamass=1.197449; 
  Double_t p1,p2,E1,E2;
  Int_t pion[2]={-1,-1}; Int_t kminus=-1; Int_t omeindx=-1;
  Int_t j=0,k=0;
  for(Int_t i=0; i<RN_nt; i++) {// particle IDs pi+:8 pi-:9 k-:12 p:14
    if(mvert[i]==1) {           // particle IDs pi+:8 pi-:9 k-:12 p:14 Ome_minus=24 X_minus=23
      if(pid[i]==8 || pid[i]==9) {pion[j]=i;j++;}
      if(pid[i]==24) omeindx=i;
    }
  } 
  if(pion[1]!=-1) cout<<"there is a second pi+!!!!"<<endl;
  // now ReCalculate the Omegac0 mass 
  p1=sqrt(pow(p_x[pion[0]],2) + pow(p_y[pion[0]],2) + pow(p_z[pion[0]],2)); //  pi+
  p2=sqrt(pow(p_x[omeindx],2) + pow(p_y[omeindx],2) + pow(p_z[omeindx],2)); // Omega- 
  E1=sqrt(pow(m[pion[0]],2) + pow(p1,2)); 
  //E2=sqrt(pow(m[omeindx],2) + pow(p2,2));  // 
  //cout<<m[piplus[0]] <<" "<<p1<<" "<<E1<<endl;
  E2=sqrt(pow(Ximass,2) + pow(p2,2));  // Replace Omega- with Xi-
  Double_t newmass= sqrt(pow(E1+E2,2) -
			 pow(p_x[pion[0]]+p_x[omeindx],2) -
			 pow(p_y[pion[0]]+p_y[omeindx],2) -
			 pow(p_z[pion[0]]+p_z[omeindx],2)
			 );
  //cout<<"newmass"<<newmass<<endl;
  return newmass;
}


Double_t s::NewMass(Int_t entry) { // written for 850
   Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
                     RN_pidt8, RN_pidt9} ;
   Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
                       RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
   Float_t mass[9] = {RN_masst1, RN_masst2, RN_masst3, RN_masst4, RN_masst5, RN_masst6, 
                      RN_masst7, RN_masst8, RN_masst9} ;
   Int_t no_of_trk = RN_nt;
   Double_t p_x[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9};
   Double_t p_y[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9};
   Double_t p_z[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9};
  
  Double_t mass_ome=1.67243; Double_t mass_xi=1.32131; Double_t mass_pi=0.139568; 
  // Double_t mass_ome=1.67245; Double_t mass_xi=1.32131; Double_t mass_pi=0.13957018; 
  
  Double_t E_pi, E_ome; //Energies
  Double_t P_pi, P_ome; //Momenta
  Int_t pi_indx=-1; //
  Int_t omeindx=-1;
  for(Int_t i=0; i<no_of_trk; i++) {// particle IDs pi+:8 pi-:9 k-:12 p:14
    if(pid[i]==8  && mvert[i]==1) pi_indx=i; //found the pion plus
    if(pid[i]==24 && mvert[i]==1) omeindx=i; 
  }

  P_pi=sqrt(pow(p_x[pi_indx],2)+ pow(p_y[pi_indx],2)+ pow(p_z[pi_indx],2)); // pion momentum
  E_pi=sqrt(pow(mass[pi_indx],2) + pow(P_pi,2)); //pion energy

  P_ome=sqrt(pow(p_x[omeindx],2) + pow(p_y[omeindx],2) + pow(p_z[omeindx],2)); // omega momentum
  //E_ome=sqrt(pow(mass[omeindx],2) + pow(P_ome,2)); //omega energy
  E_ome=sqrt(pow(mass_xi,2) + pow(P_ome,2)); // Xi- energy
  
  
  //  Double_t m1=mass_ome;  // 1 : Omega-
  Double_t m1=mass_xi;  // 1 : Replace Ome- with Xi-
  Double_t m2=mass_pi;// 2 : pion
  Double_t E1=E_ome;
  Double_t E2=E_pi;
  Double_t p1x=p_x[omeindx]; Double_t p1y=p_y[omeindx]; Double_t p1z=p_z[omeindx];
  Double_t p2x=p_x[pi_indx]; Double_t p2y=p_y[pi_indx]; Double_t p2z=p_z[pi_indx];
  
  
  Double_t newmass = sqrt (pow(m1+m2,2) + 2*(E1*E2 - m1*m2 - (p1x*p2x + p1y*p2y + p1z*p2z) ) );
  
  return newmass;
  

}
void s::Order(Double_t *p, Double_t *q) {//return p,q in ascending order
  int temp;
  if(*p>*q) {
	temp=*p;
	*p=*q;
	*q=temp;
  }
}
void s::Bubble(Double_t *a,int n) {//Bubble sorting of float array A[]
  int i,j;
  
  for (i=0; i<n; i++){
    for (j=n-1; i<j; j--){
      Order(&a[j-1], &a[j]);
    }
  }
}

Double_t s::Drmass(Int_t entry) {
  //cout<<RN_X0t1<<" "<<RN_Y0t1<<" "<<RN_Pxt1<<" "<<RN_Pyt1<<" "<<RN_Pzt1<<" "<<RN_pidt1<<" "<<RN_chit1<<" "<<RN_pmtrkt1<<endl;
  Float_t tt[8][9]={{RN_X0t1,RN_X0t2,RN_X0t3,RN_X0t4,RN_X0t5,RN_X0t6,RN_X0t7,RN_X0t8,RN_X0t9},
		    {RN_Y0t1,RN_Y0t2,RN_Y0t3,RN_Y0t4,RN_Y0t5,RN_Y0t6,RN_Y0t7,RN_Y0t8,RN_Y0t9},
		    {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9},
		    {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9},
		    {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9},
		    {RN_pidt1,RN_pidt2,RN_pidt3,RN_pidt4,RN_pidt5,RN_pidt6,RN_pidt7,RN_pidt8,RN_pidt9},
		    {RN_chit1,RN_chit2,RN_chit3,RN_chit4,RN_chit5,RN_chit6,RN_chit7,RN_chit8,RN_chit9},
		    {RN_pmtrkt1,RN_pmtrkt2,RN_pmtrkt3,RN_pmtrkt4,RN_pmtrkt5,RN_pmtrkt6,RN_pmtrkt7,RN_pmtrkt8,RN_pmtrkt9}};
		    
  

  //cout<<RN_pmtrkt1<<" "<<RN_pmtrkt2<<" "<<RN_pmtrkt3<<" "<<RN_pmtrkt4<<RN_pmtrkt5<<RN_pmtrkt6<<RN_pmtrkt7<<RN_pmtrkt8<<RN_pmtrkt9<<endl;
  //Int_t pmt[9]={RN_pmtrkt1,RN_pmtrkt2,RN_pmtrkt3,RN_pmtrkt4,RN_pmtrkt5,RN_pmtrkt6,RN_pmtrkt7,RN_pmtrkt8,RN_pmtrkt9}; //equivalence pmtrk
  Int_t pmt[9]={RN_mvert_t1,RN_mvert_t2,RN_mvert_t3,RN_mvert_t4,RN_mvert_t5,RN_mvert_t6,RN_mvert_t7,RN_mvert_t8,RN_mvert_t9}; //equivalence pmtrk

  Double_t masses[24]={0., 0.00051099, 0.00051099, 0., 0.105658389, 0.105658389,
		       0.1349743, 0.1395679, 0.1395679, 0.497671, 0.493646, 0.493646, 
		       0.93956563, 0.93827231, 0.93827231,0.497671, 0.54745, 
		       1.11563, 1.18937, 1.19255, 1.197465, 1.31485, 1.32133, 1.67243};
  // 8 pion+          0.1395679      1       2.603e-8      0
  // 9 pion-          0.1395679     -1       2.603e-8      0
  //11 kaon+          0.493646       1       1.237e-8      0
  //12 kaon-          0.493646      -1       1.237e-8      0
  //14 proton         0.93827231     1       1.0e+15       0
  //15 antiproton     0.93827231    -1       1.0e+15       0

  Double_t rmass; //local mass
  Double_t p[4]={0,0,0,0};  // local 4-vector
  Float_t kk;             //  track limits
  Double_t c[8][20];       //  track measurements, 4-vec, mass
  Float_t pl,cc,tx,ty,ee; //  local track measurements, 4-vec components
  Double_t dmdc[3];        //  mass - measurement derivative
  Double_t dpdc[4][3];     //  4-vector component - measurement derivative
  Double_t sig[3][3][20];  //  track measurement error matrix
  Float_t s;              //  momentum scale factor
  Double_t t;              //  sum variable
  Int_t md;                 //  particle mass id
  Int_t itgt;               //  upstream target
  Int_t itrk,jtrk,i,j,k,l,n;  //  dummy indices
  //  angular resolution parameters
  Double_t stx_par[3][6]={{5.307e-6, 7.268e-6, 8.499e-6, 9.868e-6, 12.968e-6, 16.069e-6},
			  {4763.e-6, 3264.e-6, 2828.e-6, 2315.e-6, 1658.e-6, 1176.e-6},
			  {-166.e-6, -231.e-6, -275.e-6, -343.e-6, -530.e-6, -872.e-6}};
  //  momentum resolution parameters
  Double_t sp_par[2][2]={{96.e-6,  34.e-6},  //      ! M1 tracks
			 {11.4,  23.0}};     //      ! M2 tracks
  //-------------------------------------------------------------------------------------------------  
  Double_t drmass=0.; //default value
  for (int ii=0;ii<4;ii++) {p[ii]=0;} //initialize p 
  
  for( int itrk = 0; itrk<RN_nt; itrk++) { //   each track
    jtrk    = pmt[itrk] ;  //   track index
    if (jtrk==1) {                           //    which is in vtx
      //      md        = (pmt[itrk] & 255);        //    particle number (ID)
      md        = tt[5][itrk];              //    particle number (ID)
      rmass     = masses[8];                //    default mass - pi-
      if (md>0) rmass = masses[md-1];       //    get mass ,  ome-=24, but it is index starts from 0 !!!
                                            //   so masses[md-1]
      c[0][itrk] = tt[2][itrk]/tt[4][itrk]; //    track measurements
      c[1][itrk] = tt[3][itrk]/tt[4][itrk];
      c[2][itrk] = 1./sqrt( pow(tt[2][itrk],2) + pow(tt[3][itrk],2) + pow(tt[4][itrk],2)); 
      
      for(int ii=0;ii<20;ii++){sig[1][1][ii]=0;} // initia;ize sig[1][1][itrk]
      
      itgt            = TMath::Max(1,TMath::Min(6,int((RN_Zp+7.5)/1.5)));

      sig[0][0][itrk]   =pow((stx_par[0][itgt] 
			      + stx_par[1][itgt]*c[2][itrk] 
			      + stx_par[2][itgt]*pow(c[2][itrk],2)),2) ;
      sig[0][1][itrk]   = 0.0;
      sig[1][1][itrk]   = sig[0][0][itrk] ;
      
      j                 = 0;
      if (c[2][itrk] < 1./15.) j=1;
      sig[0][2][itrk]   = 0.0;
      sig[1][2][itrk]   = 0.0; //curvature errors
      sig[2][2][itrk]   = pow(sp_par[0][j],2) *  (1.0+ pow((sp_par[1][j]*c[2][itrk]),2)) ;

      sig[1][0][itrk]   = sig[0][1][itrk] ;     // ! diagonalize matrix
      sig[2][0][itrk]   = sig[0][2][itrk] ;
      sig[2][1][itrk]   = sig[1][2][itrk] ;
      
      c[3][itrk] = sqrt( pow(rmass,2) + pow(tt[2][itrk],2) + 
			 pow(tt[3][itrk],2) + pow(tt[4][itrk],2)); //! E
      c[4][itrk] = tt[2][itrk] ;                                   //  Px
      c[5][itrk] = tt[3][itrk] ;                                   //  Py
      c[6][itrk] = tt[4][itrk] ;                                   //  Pz
      c[7][itrk] = rmass ;                                         //  mass
    
      for(int n=0; n<4;n++) {                  //  4-vector components
	p[n]   = p[n] + c[3+n][itrk];           //  sum 4 vector 
      }                                         //  4-vector components
      // recon_mass.F version 
      //p[3] = abs(1./sqrt(1+ pow(c[0][itrk],2))/c[2][itrk]);            // z momentum component
      //p[2] = p[3]*c[1][itrk];                                          // y momentum component
      //p[1] = p[3]*c[0][itrk] ;                                         // x momentum component
      //p[0] = sqrt(pow(rmass,2)+pow(p[1],2)+pow(p[2],2)+pow(p[3],2));   // energy
    }     //  which is in vtx
  }       //  each track
  
  Double_t rmass     = sqrt(TMath::Max(0., pow(p[0],2) - pow(p[1],2) -pow(p[2],2) - pow(p[3],2))) ;
  
     t = 0.;                                //  sum variable
 
   for(itrk=0;itrk<RN_nt;itrk++) {             //  each vertex track   :     
    jtrk    = pmt[itrk] ;
    
    if(jtrk==1) {
      tx          = c[0][itrk];                 //   track x slope
      ty          = c[1][itrk];                 //   track y slope
      cc          = c[2][itrk];                 //   track curvature (1/sqrt(Px^2_Pz^2)
      ee          = c[3][itrk];                 //   track energy
      pl          = c[6][itrk];                 //   track z momentum component
      //   full derivatives 
      dpdc[0][0]   = -tx *  pow(ty,2)  / (ee * pow(((1+pow(tx,2))*cc),2)); //! dE /dTx
      dpdc[1][0]   =  pl * (1-pow(tx*pl*cc,2)) ;             // dPx/dTx
      dpdc[2][0]   = - pow(pl,3) * pow(cc,2) * tx *ty ;      // dPy/dTx
      dpdc[3][0]   = - pow(pl,3) * pow(cc,2) * tx;           // dPz/dTx
      
      dpdc[0][1]   =  ty/(ee*(1+ pow(tx,2)) * pow(cc,2)) ;   // dE /dTy
      dpdc[1][1]   =  0.0;                                   // dPx/dTy
      dpdc[2][1]   =  pl;                                    // dPy/dTy
      dpdc[3][1]   =  0.0;                                   // dPz/dTy
      
      dpdc[0][2]   = -(1+ pow(ty,2)/(1+ pow(tx,2)))/(ee * pow(cc,3)); // dE /dc
      dpdc[1][2]   = -pl*tx/cc ;                                      // dPx/dc
      dpdc[2][2]   = -pl*ty/cc;                                       // dPy/dc
      dpdc[3][2]   = -pl/cc ;                                         // dPz/dc

      //Double_t dmdc[3]={0.,0.,0.}; //initialize dmdc[i]
      
      for( i=0;i<3;i++) {                         // each variable
	dmdc[i] = p[0] * dpdc[0][i];              // do dot products
	for( l=1;l<4;l++) {                 
	  dmdc[i] = dmdc[i] - p[l] * dpdc[l][i];          
	}
	dmdc[i]   = dmdc[i]/rmass;
      }                                        // each variable
      
      for( i=0;i<3;i++) {                        // propagate error
	for( j=0;j<3;j++) {
	  t = t + dmdc[i] * dmdc[j] * sig[i][j][itrk];
	}
      }
      // cout<<"t "<<t<<endl;
      
    } // if(jtrk==1)
   }    // each track
   //if(t<1.e-300) cout<< "t "<<t<<"sig "<<sig[2][2][itrk]<<"drmass "<<drmass<<endl;
   t= TMath::Max(t,1.e-8);
   drmass=sqrt(t);
   //if(drmass>0.05) cout<<drmass<<endl;
   return drmass;
    // 
    
}
 
Int_t s::FindParticle(Int_t entry, const char *part_name) { 
  //return asked particles track index...
  
  Float_t m[9] = {RN_masst1,RN_masst2,RN_masst3,RN_masst4,RN_masst5,RN_masst6,RN_masst7,
		  RN_masst8,RN_masst9} ;
  Float_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Int_t pid1[9] = {RN_pid1t1, RN_pid1t2, RN_pid1t3, RN_pid1t4, RN_pid1t5, RN_pid1t6, RN_pid1t7,
		   RN_pid1t8, RN_pid1t9} ;
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
  Float_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
  Float_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ;
  Float_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9} ;
  Int_t i,j=0; // particle IDs pi+:8 pi-:9 k-:12 p:14 Ome-:24 X-:23  sigma +:19, sigma -:21
  // particle IDs pi+:8 pi-:9 k-:12 p:14 Ome_minus=24 X_minus=23
 if (part_name=="k-")  part_id=12;      // k- =12, 
 if (part_name=="pi+") part_id=8;      //pi+ =8,
 if (part_name=="pi-") part_id=9;     // pi- =9, 
 if (part_name=="p")   part_id=14;   //p  =14
 if (part_name=="ome-") part_id=24; // ome-  =24

  for(Int_t i=0; i<RN_nt ; i++) {
    if(RN_pid[i]==part_id) return i;
  }
}


Double_t s::DaughterMom(Int_t entry) {
  Float_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;

  Float_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
  Float_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ;
  Float_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9} ;

  for(Int_t i=0; i<RN_nt ; i++) {
    if(mvert[i]==2) {Double_t momentum=Trk_p(entry,i); return momentum;}
  }
}

void s::hex2bin(Int_t n) {

  //while (cin >> n) {
  //    cout << "decimal: " << n << endl;
    // print binary with leading zeros
    cout << "binary : ";
    for (int i=31; i>=0; i--) {
      int bit = ((n >> i) & 1);
	cout << bit;
    }
    cout << endl;
    //  }//end loop
} //void

void s::GenerateEmbedFile(Int_t entry) {
  Int_t mvert[9] = {RN_mvert_t1, RN_mvert_t2, RN_mvert_t3, RN_mvert_t4, RN_mvert_t5, 
		      RN_mvert_t6, RN_mvert_t7, RN_mvert_t8, RN_mvert_t9} ;
  Double_t px[9] = {RN_Pxt1,RN_Pxt2,RN_Pxt3,RN_Pxt4,RN_Pxt5,RN_Pxt6,RN_Pxt7,RN_Pxt8,RN_Pxt9} ;
  Double_t py[9] = {RN_Pyt1,RN_Pyt2,RN_Pyt3,RN_Pyt4,RN_Pyt5,RN_Pyt6,RN_Pyt7,RN_Pyt8,RN_Pyt9} ;
  Double_t pz[9] = {RN_Pzt1,RN_Pzt2,RN_Pzt3,RN_Pzt4,RN_Pzt5,RN_Pzt6,RN_Pzt7,RN_Pzt8,RN_Pzt9} ;
  Double_t m[9] = {RN_masst1,RN_masst2,RN_masst3,RN_masst4,RN_masst5,RN_masst6,RN_masst7,
		  RN_masst8,RN_masst9} ;
  Double_t pid[9] = {RN_pidt1, RN_pidt2, RN_pidt3, RN_pidt4, RN_pidt5, RN_pidt6, RN_pidt7, 
		    RN_pidt8, RN_pidt9} ;
  Double_t x[9] = {RN_X0t1,RN_X0t2,RN_X0t3,RN_X0t4,RN_X0t5,RN_X0t6,RN_X0t7,RN_X0t8,RN_X0t9} ;
  Double_t y[9] = {RN_Y0t1,RN_Y0t2,RN_Y0t3,RN_Y0t4,RN_Y0t5,RN_Y0t6,RN_Y0t7,RN_Y0t8,RN_Y0t9} ;
  Int_t pi_pl[3]={-1,-1,-1}; Int_t pi_mi[3]={-1,-1,-1}; Int_t ka=-1;Int_t ome=-1;
  Int_t j=0; Int_t k=0; 
  for(Int_t i=0; i<RN_nt ; i++) {
    if (pid[i]==8 && mvert[i]==1) {pi_pl[j]=i; j++;}
    else 
      if (pid[i]==9 && mvert[i]==1) {pi_mi[k]=i; k++;}
      else 
	if (mvert[i]==2) {ka=i;}
	else 
	  if (pid[i]==24 && mvert[i]==1) {ome=i; }
  }
  // track momenta
  Double_t tk1p=sqrt(px[ome]*px[ome] + py[ome]*py[ome] + pz[ome]*pz[ome] );
  Double_t tk2p=sqrt(px[pi_pl[0]]*px[pi_pl[0]] + py[pi_pl[0]]*py[pi_pl[0]] + pz[pi_pl[0]]*pz[pi_pl[0]] );
  Double_t tk3p=sqrt(px[pi_mi[0]]*px[pi_mi[0]] + py[pi_mi[0]]*py[pi_mi[0]] + pz[pi_mi[0]]*pz[pi_mi[0]] );
  Double_t tk4p=sqrt(px[pi_pl[1]]*px[pi_pl[1]] + py[pi_pl[1]]*py[pi_pl[1]] + pz[pi_pl[1]]*pz[pi_pl[1]] );
  Double_t tk5p=sqrt(px[ka]*px[ka] + py[ka]*py[ka] + pz[ka]*pz[ka]);
  // track x,y,z coordinates
  Double_t m1=1.672; Double_t x1=x[ome];       Double_t y1=y[ome];       Double_t z1=RN_Zs; 
  Double_t m2=0.139; Double_t x2=x[pi_pl[0]];  Double_t y2=y[pi_pl[0]];  Double_t z2=RN_Zs; 
  Double_t m3=0.140; Double_t x3=x[pi_mi[0]];  Double_t y3=y[pi_mi[0]];  Double_t z3=RN_Zs; 
  Double_t m4=0.140; Double_t x4=x[pi_pl[1]];  Double_t y4=y[pi_pl[1]];  Double_t z4=RN_Zs; 
  Double_t m5=0.494; Double_t x5=x[ka];        Double_t y5=y[ka];        Double_t z5=RN_Zs2; 
  // track directional cosines  
  Double_t cos1x=px[ome]/tk1p ;       Double_t cos1y=py[ome]/tk1p ;       Double_t cos1z=pz[ome]/tk1p ;
  Double_t cos2x=px[pi_pl[0]]/tk2p ;  Double_t cos2y=py[pi_pl[0]]/tk2p ;  Double_t cos2z=pz[pi_pl[0]]/tk2p ; 
  Double_t cos3x=px[pi_mi[0]]/tk3p ;  Double_t cos3y=py[pi_mi[0]]/tk3p ;  Double_t cos3z=pz[pi_mi[0]]/tk3p ; 
  Double_t cos4x=px[pi_pl[1]]/tk4p ;  Double_t cos4y=py[pi_pl[1]]/tk4p ;  Double_t cos4z=pz[pi_pl[1]]/tk4p ; 
  Double_t cos5x=px[ka]/tk5p ;        Double_t cos5y=py[ka]/tk5p ;        Double_t cos5z=pz[ka]/tk5p ; 
  // write into file
  fprintf (pFile, "%7.4f %7.4f %7.4f %8.6f %8.6f %8.6f %8.6f %5.3f\n ",x1,y1,z1,cos1x,cos1y,cos1z,tk1p,m1);
  fprintf (pFile, "%7.4f %7.4f %7.4f %8.6f %8.6f %8.6f %8.6f %5.3f\n ",x2,m2,z2,cos2x,cos2y,cos2z,tk2p,m2);
  fprintf (pFile, "%7.4f %7.4f %7.4f %8.6f %8.6f %8.6f %8.6f %5.3f\n ",x3,y3,z3,cos3x,cos3y,cos3z,tk3p,m3);
  fprintf (pFile, "%7.4f %7.4f %7.4f %8.6f %8.6f %8.6f %8.6f %5.3f\n ",x4,y4,z4,cos4x,cos4y,cos4z,tk4p,m4);
  fprintf (pFile, "%7.4f %7.4f %7.4f %8.6f %8.6f %8.6f %8.6f %5.3f\n ",x5,y5,z5,cos5x,cos5y,cos5z,tk5p,m5);
  }

void s::NoteEvent(Int_t entry) {
  fprintf(pFile, "------------------------------------------------------------------------- \n");
  fprintf(pFile, "run# %8d,  event# %12d ID# %6d  mass %8.4f\n",int(RN_run), int(RN_event),RN_ID,RN_mass);
  fprintf(pFile, "kink  usec  sip  sin  ksi  ome  oxi \n");
  Bool_t kink=false; if((RN_status_rv0&kink_s_v0)!=0) kink=true;
  Bool_t usec=false; if((RN_status_rv0&usec_s_v0)!=0) usec=true;
  Bool_t sip=false; if((RN_status_rv0&sip_s_v0)!=0)  sip=true;
  Bool_t sin=false; if((RN_status_rv0&sin_s_v0)!=0)  sin=true;
  Bool_t ksi=false; if((RN_status_rv0&ksi_s_v0)!=0)  ksi=true;
  Bool_t ome=false; if((RN_status_rv0&ome_s_v0)!=0)  ome=true;
  Bool_t oxi=false; if((RN_status_rv0&oxi_s_v0)!=0)  oxi=true;
  
  fprintf(pFile, "%4d %4d %4d %4d %4d %4d %4d\n",kink,usec,sip,sin,ksi,ome,oxi);
  fprintf(pFile, "------------------------------------------------------------------------- \n");
 
 }

void s::DumpEvent(Int_t entry) {
  //cout<<"----------------------------------------------------------------"<<endl;
  fprintf(pFile, "run# %8d,  event# %12d ID# %6d # of sec vtx: %d\n",int(RN_run), int(RN_event),RN_ID,RN_ns);
  fprintf(pFile, "            x       sx        y       sy       z       sz        chi2       nf\n")  ;
  fprintf(pFile, "pvtx   %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.2f %8d \n",RN_Xp,RN_Sxp,RN_Yp,RN_Syp,RN_Zp,RN_Szp,RN_Pchi,RN_Pnf);
  fprintf(pFile, "svtx   %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.2f \n",RN_Xs,RN_Sxs,RN_Ys,RN_Sys,RN_Zs,RN_Szs,RN_Schi);
  fprintf(pFile, "svtx2  %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.2f \n",RN_Xs2,RN_Sxs2,RN_Ys2,RN_Sys2,RN_Zs2,RN_Szs2,RN_Schi2);
  fprintf(pFile, "            x       px        y       py       pz      chi2   mvert    pid    pmtrk\n")  ;
  fprintf(pFile, "beam   %8.4f %8.4f %8.4f %8.4f %8.4f\n",RN_x0b,RN_Pxb,RN_y0b,RN_Pyb,RN_Pzb);
  fprintf(pFile, "tk1    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t1,RN_Pxt1,RN_Y0t1,RN_Pyt1,RN_Pzt1,RN_chit1,RN_mvert_t1,RN_pidt1,RN_pmtrkt1);
  fprintf(pFile, "tk2    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t2,RN_Pxt2,RN_Y0t2,RN_Pyt2,RN_Pzt2,RN_chit2,RN_mvert_t2,RN_pidt2,RN_pmtrkt2);
  fprintf(pFile, "tk3    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t3,RN_Pxt3,RN_Y0t3,RN_Pyt3,RN_Pzt3,RN_chit3,RN_mvert_t3,RN_pidt3,RN_pmtrkt3);
  fprintf(pFile, "tk4    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t4,RN_Pxt4,RN_Y0t4,RN_Pyt4,RN_Pzt4,RN_chit4,RN_mvert_t4,RN_pidt4,RN_pmtrkt4);

  if(RN_nt==5){
    fprintf(pFile, "tk5    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t5,RN_Pxt5,RN_Y0t5,RN_Pyt5,RN_Pzt5,RN_chit5,RN_mvert_t5,RN_pidt5,RN_pmtrkt5);
  }
  if(RN_nt==6){
    fprintf(pFile, "tk5    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t5,RN_Pxt5,RN_Y0t5,RN_Pyt5,RN_Pzt5,RN_chit5,RN_mvert_t5,RN_pidt5,RN_pmtrkt5);
    fprintf(pFile, "tk6    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t6,RN_Pxt6,RN_Y0t6,RN_Pyt6,RN_Pzt6,RN_chit6,RN_mvert_t6,RN_pidt6,RN_pmtrkt6);
  }
   if(RN_nt==7){
    fprintf(pFile, "tk5    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t5,RN_Pxt5,RN_Y0t5,RN_Pyt5,RN_Pzt5,RN_chit5,RN_mvert_t5,RN_pidt5,RN_pmtrkt5);
    fprintf(pFile, "tk6    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t6,RN_Pxt6,RN_Y0t6,RN_Pyt6,RN_Pzt6,RN_chit6,RN_mvert_t6,RN_pidt6,RN_pmtrkt6);
    fprintf(pFile, "tk7    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t7,RN_Pxt7,RN_Y0t7,RN_Pyt7,RN_Pzt7,RN_chit7,RN_mvert_t7,RN_pidt7,RN_pmtrkt7);
   }
   if(RN_nt==8){
     fprintf(pFile, "tk5    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t5,RN_Pxt5,RN_Y0t5,RN_Pyt5,RN_Pzt5,RN_chit5,RN_mvert_t5,RN_pidt5,RN_pmtrkt5);
     fprintf(pFile, "tk6    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t6,RN_Pxt6,RN_Y0t6,RN_Pyt6,RN_Pzt6,RN_chit6,RN_mvert_t6,RN_pidt6,RN_pmtrkt6);
     fprintf(pFile, "tk7    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t7,RN_Pxt7,RN_Y0t7,RN_Pyt7,RN_Pzt7,RN_chit7,RN_mvert_t7,RN_pidt7,RN_pmtrkt7);
     fprintf(pFile, "tk8    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t8,RN_Pxt8,RN_Y0t8,RN_Pyt8,RN_Pzt8,RN_chit8,RN_mvert_t8,RN_pidt8,RN_pmtrkt8);
   }
   if(RN_nt==9){
     fprintf(pFile, "tk5    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t5,RN_Pxt5,RN_Y0t5,RN_Pyt5,RN_Pzt5,RN_chit5,RN_mvert_t5,RN_pidt5,RN_pmtrkt5);
     fprintf(pFile, "tk6    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t6,RN_Pxt6,RN_Y0t6,RN_Pyt6,RN_Pzt6,RN_chit6,RN_mvert_t6,RN_pidt6,RN_pmtrkt6);
     fprintf(pFile, "tk7    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t7,RN_Pxt7,RN_Y0t7,RN_Pyt7,RN_Pzt7,RN_chit7,RN_mvert_t7,RN_pidt7,RN_pmtrkt7);
     fprintf(pFile, "tk8    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t8,RN_Pxt8,RN_Y0t8,RN_Pyt8,RN_Pzt8,RN_chit8,RN_mvert_t8,RN_pidt8,RN_pmtrkt8);
     fprintf(pFile, "tk9    %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8d %8x\n",RN_X0t9,RN_Pxt9,RN_Y0t9,RN_Pyt9,RN_Pzt9,RN_chit9,RN_mvert_t9,RN_pidt9,RN_pmtrkt9);
   }
   fprintf(pFile, "     mass     px        py        pz         e         L       sigma\n");
   fprintf(pFile, " %8.4f  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f\n",RN_mass,RN_Pxs, RN_Pys,RN_Pzs,RN_Es,RN_L,RN_S);
   fprintf(pFile, "     mass2   px2        py2       pz2        e2\n");
   fprintf(pFile, " %8.4f  %8.4f  %8.4f  %8.4f  %8.4f\n",RN_mass2,RN_Pxs2, RN_Pys2,RN_Pzs2,RN_Es2);

   fprintf(pFile, "          ********** Svkink Objects *********\n"); 
   fprintf(pFile, "  chi_rv02  dca_rv02 tet_rv02\n"); 
   fprintf(pFile, "%8.4f  %8.4f  %8.4f\n",RN_chi_rv02,RN_dca_rv02,RN_tet_rv02); 
   fprintf(pFile, "  \n");
   fprintf(pFile, "  pd_ph   chi_rv01 dca_rv01 pidtrk_rv01\n"); 
   fprintf(pFile, "%8.4f %8.4f %8.4f %8d\n",RN_pd_ph,RN_chi_rv01,RN_dca_rv01,RN_pidtrk_rv01); 
   fprintf(pFile, "  \n");
  
   fprintf(pFile, "  pd_ph_b chi_rv0b dca_rv0b z_rv0b     vid_rv01b vid_rv02b status_rv0\n");
   fprintf(pFile, "%8.4f %8.4f %8.4f %9.4f %10d %10d %10x\n",RN_pd_ph_b,RN_chi_rv0b,RN_dca_rv0b,RN_z_rv0b,RN_vid_rv01b,RN_vid_rv02b,RN_status_rv0);
   
   fprintf(pFile, "kink  usec  sip  sin  ksi  ome  oxi \n");
   Bool_t kink=false; if((RN_status_rv0&kink_s_v0)!=0) kink=true;
   Bool_t usec=false; if((RN_status_rv0&usec_s_v0)!=0) usec=true;
   Bool_t sip=false; if((RN_status_rv0&sip_s_v0)!=0)  sip=true;
   Bool_t sin=false; if((RN_status_rv0&sin_s_v0)!=0)  sin=true;
   Bool_t ksi=false; if((RN_status_rv0&ksi_s_v0)!=0)  ksi=true;
   Bool_t ome=false; if((RN_status_rv0&ome_s_v0)!=0)  ome=true;
   Bool_t oxi=false; if((RN_status_rv0&oxi_s_v0)!=0)  oxi=true;

   fprintf(pFile, "%4d %4d %4d %4d %4d %4d %4d\n",kink,usec,sip,sin,ksi,ome,oxi);
   fprintf(pFile, "------------------------------------------------------------------------- \n");
 
   
}

#endif // #ifdef s_cxx