{

AutoTree->Process("Selector.C")
}

#define Selector_cxx
// The class definition in Selector.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.


// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// root> T->Process("Selector.C")
// root> T->Process("Selector.C","some options")
// root> T->Process("Selector.C+")
//


#include "Selector.h"
#include <TH2.h>
#include <TStyle.h>

//================================SETTINGS ========
const int chanmax=8000;
const int bins2d=800;

//const long eventsmax=2e+10;
const long eventsmax=1e+6;

//calibration-------------EXL   19 - 1st is zero

double a[19]={0.,
0.668983,0.519247,0.5649,
0.600488,0.579928,0.598377,
0.540216,0.488096,0.610848,
0.560869, 0.503247,0.464832,
0.848081,0.792648,0.799,	
0.426236,0.485055,0.492694};

double b[19]={0.,
 -71.1150, 132.36, 71.2,
 -58.9582, 72.8547, -56.2703,
 91.2014, 95.1007, 9.26351,
 -2.42136, 86.8976,  105.457,
 45.8667, 68.9514,2.84729,
 18.4952, 18.7380, 116.032};

//thresholds EXL
double tr[19]={0,
140.,140.,140.,
140.,140.,140.,
140.,140.,140.,
140.,140.,140.,
140.,140.,140.,
140.,140.,140.
};


// Experimental factor R1ex from MTX: includes gain G1
//  the theoretical factor R1 should be used in correction

// higher positive slope line - like exl1,3,5 signals
double R1ex[9]={
 0.14,0.18,0.15,
 0.16,0.18,0.16, 
 0.15,0.42,0.17
};
// lower positive slope line - like exl2,4,6 signal
double R2ex[9]={
 0.16,0.15,0.14,
 0.15,0.15,0.16, 
 0.16,0.15,0.16
};

//--- we compute gain ratio later
double G1comp[9];


// arbitrarym good around 1.0
// we recalculate later
double G2gain[9]={
  2616., 2368., 2639.,
  2717.,  2728., 4308., 
  1822.,  3168., 2796.
};
double G2gain_factor=3000.;

//=============================SETTTING   END _-------



void Selector::Begin(TTree * /*tree*/)
{
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).

  char ch1[100];
  char ch2[100];
   TString option = GetOption();
   printf("beginning the process\n%s","");

   //  HISTO DEFINITIONS-------------------

   for (int i=0;i<18;i++){
     //  SINGLES
     sprintf(ch1, "exl%02d", i+1 );
     sprintf(ch2, "%s   SINGLES", ch1 );
     exl[i]=new TH1F(ch1,ch2,chanmax,0,chanmax);
   }
   for (int i=0;i<18;i++){
     //  SINGLES CALIB
     sprintf(ch1, "exlcal%02d", i+1 );
     sprintf(ch2, "%s   SINGLES CALIB", ch1 );
     exlcal[i]=new TH1F(ch1,ch2,chanmax,0,chanmax);
   }

   exl_singlesum=new TH1F("exlcalSUM","TOTAL SUM of SINGLES / useful4 online",chanmax,0,chanmax);   
   exl_pairsum=new TH1F("pairSUM","TOTAL SUM of summed pairs/ wrong concept",chanmax,0,chanmax);

   complete_pairs=new TH1F("complpairs","pair hits in the matrix",18,0,18);
   incomplete_pairs=new TH1F("cincomppairs","single hits in pairs in the matrix",18,0,18);

   // BIDIM--------------------------------------


   for (int i=0;i<9;i++){
     //  PAIRS
     int j1=i*2+1;
     int j2=i*2+2;
     sprintf(ch1, "pair%02d%02d",  j1, j2);
     sprintf(ch2, "%s   PAIRS", ch1 );
     expair[i]=new TH2F(ch1,ch2,bins2d,0,chanmax,bins2d,0,chanmax);
   }
   for (int i=0;i<9;i++){
     //  PAIRS  CALIB
     int j1=i*2+1;
     int j2=i*2+2;
     sprintf(ch1, "paircal%02d%02d",  j1, j2);
     sprintf(ch2, "%s   PAIRS CALIB", ch1 );
     expaircal[i]=new TH2F(ch1,ch2,bins2d,0,chanmax,bins2d,0,chanmax);
   }
   for (int i=0;i<9;i++){
     //  SUMPAIRS SINGLE 
     int j1=i*2+1;
     int j2=i*2+2;
     sprintf(ch1, "paircalsum%02d%02d",  j1, j2);
     sprintf(ch2, "%s   SINGLES = PAIRS CALIBrated AND SUMMED, energy restored", ch1 );
     exlsumpair[i]=new TH1F(ch1,ch2,chanmax,0,chanmax);
   }




   for (int i=0;i<6;i++){
     //  HNEIGH
     int j1=i*2+2;
     if (j1>=6){ j1=j1+2;}
     if (j1>=12){ j1=j1+2;}
     sprintf(ch1, "hneig%02d%02d", j1,j1+1 );
     sprintf(ch2, "%s   HNEIGHBOR", ch1 );
     exhneig[i]=new TH2F(ch1,ch2,bins2d,0,chanmax,bins2d,0,chanmax);
   }
   for (int i=0;i<12;i++){
     //  VNEIGH
     int j1=i+1;
     sprintf(ch1, "vneigh%02d%02d", j1,j1+6 );
     sprintf(ch2, "%s   VNEIGHBOR", ch1 );
     exvneig[i]=new TH2F(ch1,ch2,bins2d,0,chanmax,bins2d,0,chanmax);
   }


   //compensation for GAIN1 - odd crystals
   for (int i=0;i<9;i++){
     // 2/1 bad sums
     G1comp[i]=TMath::Sqrt( R2ex[i]/R1ex[i] )*1.;  // G2==1
     //     G1comp[i]=TMath::Sqrt( R1ex[i]/R2ex[i] )*1.;  // G2==1 wrong

     printf("%d   %f\n",i,G1comp[i] );
   }

} //==========BEGIN END-------------------------------------------











void Selector::SlaveBegin(TTree * /*tree*/)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();

}







// invert matrix
//  for every event from PAIR matrix :  restore the energy
//-------------------------
//    
//
Bool_t Selector::inverse_mtx(  int i9 ){
  // i take raw channel
  int ii=2*i9+1;
  int jj=2*i9+2;
  double e1=EXL->GetEnergyCrystalNumber( ii );
  double e2=EXL->GetEnergyCrystalNumber( jj );
  if (e1<=tr[ii]) return kFALSE;
  if (e2<=tr[jj]) return kFALSE;

  // we try to substract the offset first
  //   a*k+b     c*(k-d) =>  b/a=d
  e1=e1-b[ii]/a[ii];
  e2=e2-b[jj]/a[jj];

  // determine Gain1 from alpha=beta
  double G2=G2gain[i9]/G2gain_factor;
  double G1=TMath::Sqrt( R2ex[i9]/R1ex[i9] )* G2;
  double R= R2ex[i9]*G2/G1; // get to the R1==R2 (alpha==beta)
  double alpha=R/(R+1.0);// from Zamoram /easy
  double beta=alpha;     // both axes have the same leak
  double fact=G1*(1-alpha)*G2*(1-beta)-G1*G2*alpha*beta;
  fact=1.0/fact;
  //matrix"
  double N1=fact*(G2*(1-beta)*e1 - G1*alpha *e2);
  double N2=fact*(G1*(1-alpha)*e2 - G2*beta *e1);

  expaircal[i9]->Fill( N1,N2 );

  return kTRUE;
}





Bool_t Selector::is_in_pair(  int ii, int jj ){
  double e1=EXL->GetEnergyCrystalNumber( ii );
  double e2=EXL->GetEnergyCrystalNumber( jj );
  if (e1<=tr[ii]) return kFALSE;
  if (e2<=tr[jj]) return kFALSE;

  // e1=a[ii]+b[ii];
  // e2=a[jj]+b[jj];

  return kTRUE;
}



Bool_t Selector::Process(Long64_t entry)
{
   // The Process() function is called for each entry in the tree (or possibly
   // keyed object in the case of PROOF) to be processed. The entry argument
   // specifies which entry in the currently loaded tree is to be processed.
   // When processing keyed objects with PROOF, the object is already loaded
   // and is available via the fObject pointer.
   //
   // This function should contain the \"body\" of the analysis. It can contain
   // simple or elaborate selection criteria, run algorithms on the data
   // of the event and typically fill histograms.
   //
   // The processing can be stopped by calling Abort().
   //
   // Use fStatus to set the return value of TTree::Process().
   //
   // The return value is currently not used.


  double e,e1,e2, f,g;


   fReader.SetEntry(entry);
   count=count+1;	
  if (count> eventsmax){
    if (count==eventsmax+1)printf("Free ride now%s\n","");
    //Abort("just abort");