#include <fstream>
#include <sstream>
#include <string>
#include <iostream>
#include <TTree.h>
#include <TClassEdit.h>
#include <TObject.h>
#include <TList.h>
#include <TSystem.h>
#include <TSystemDirectory.h>
#include <TString.h>
#include <TFile.h>
#include <TClass.h>
#include <TCanvas.h>
#include <TGraphErrors.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2D.h>
#include <TLatex.h>
#include <TPad.h>
#include <TPaveLabel.h>
#include <TChain.h>
#include <TSpline.h>
#include <TGraph.h>
#include <TMultiGraph.h>
#include <TLine.h>
#include <TStyle.h>
#include <TF1.h>
#include <TH1F.h>
#include <TTree.h>
#include <TCanvas.h>
#include <PamLevel2.h>
#include <RunInfo.h>
#include <OrbitalInfo.h>
#include <TStopwatch.h>
#include <iomanip>

#include <CaloNuclei.h>
#include <ToFNuclei.h>
//#include <ToFNaNuclei.h>

#include </gpfs/wizard/flight/analysis/menn/devel/10th/calo/flight/CaloNuclei/inc/CaloNuclei.h>
//#include </gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFNuclei/inc/ToFNuclei.h>
//#include </gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFNaNuclei/inc/ToFNaNuclei.h>
#include </gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFPatch_v5/inc/ToFPatch.h>

#include </gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFdEdx_patch/inc/ToFdEdx_patch.h>

using namespace std;
//
//


void monitor(TString dir, TString list, TString outputname, Int_t iyear, TString ioptions, Int_t nmax=500000000){

const char* alg = "STD"; // or "EXT" "STD" "EXTF" "NUC" "NUCEXT" "NUCEXTF"
const char* alg1 = "NUC"; // or "EXT" "STD" "EXTF" "NUC" "NUCEXT" "NUCEXTF"

const char* tri_or_bi = "TRI";  // gives tri-scale calibartion, "BI"  gives bi-scale
 
char htitle[50],info[50];
char *s = new char[1];

Double_t  tmin,tmax,tint;
ULong64_t tmin1,tmax1;


// 2006
if (iyear==2006) {
tmin = 1152250000 ;
tmax = 1167606000 ;
}

// 2007
if (iyear==2007) {
tmin = 1167606000 ;
tmax = 1199142000 ;
}

// 2008
if (iyear==2008) {
tmin = 1199142000;
tmax = 1230764400;
}

// 2009
if (iyear==2009) {
tmin = 1230764400;
tmax = 1262300400;
}

// 2010
if (iyear==2010) {
tmin = 1262300400;
tmax = 1293836400;
}

// 2011
if (iyear==2011) {
tmin = 1293836400;
tmax = 1325372400;
}

// 2012
if (iyear==2012) {
tmin = 1325372400;
tmax = 1356994800;
}

// 2013
if (iyear==2013) {
tmin = 1356994800;
tmax = 1388530800;
}

// 2014
if (iyear==2014) {
tmin = 1388530800;
tmax = 1420066800;
}

// 2015
if (iyear==2015) {
tmin = 1420066800;
tmax = 1451602800;
}

// Andrey definition:
if (iyear==0) {
//   1.1.2006:  1136070000
// 31.12.2015:  1451602800
tmin = 1136070000;
tmax = 1451602800;
                   }


//tmin = tmin - 10;
tmin1=tmin;

//tmax = tmax + 10;
tmax1=tmax;


cout<<"year "<<iyear<<"   tmin "<<tmin1<<"    tmax "<<tmax1<<endl;


Int_t  icount = 0;

//------------------------------------------------------------------

  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(000);
  gStyle->SetOptFit(111);
  gStyle->SetFuncColor( kRed );


        TH2F *sc_adc_he[48];

        TH2F *sc_dedx_p1[48];
        TH2F *sc_dedx_he1[48];
        TH2F *sc_dedx_c1[48];
        TH2F *sc_dedx_p2[48];
        TH2F *sc_dedx_he2[48];
        TH2F *sc_dedx_c2[48];

        TH1F *dedx_he_trk = new TH1F("dedx_he_trk","",100.,0.,15.);

        TH1F *adc_he[48];
        for(int i=0;i<48;i++) {
        sprintf(htitle,"adc_%d",(i+1)); 
        adc_he[i] = new TH1F(htitle,htitle,100,0.,1500.);
                              }


        TH1F *dedx_p1[48];
        for(int i=0;i<48;i++) {
        sprintf(htitle,"dedx_p1_%d",(i+1)); 
        dedx_p1[i] = new TH1F(htitle,htitle,100,0.,5.);
                              }

        TH1F *dedx_he1[48];
        for(int i=0;i<48;i++) {
        sprintf(htitle,"dedx_he1_%d",(i+1)); 
        dedx_he1[i] = new TH1F(htitle,htitle,100,0.,15.);
                              }

        TH1F *dedx_c1[48];
        for(int i=0;i<48;i++) {
        sprintf(htitle,"dedx_c1_%d",(i+1)); 
        dedx_c1[i] = new TH1F(htitle,htitle,100,10.,100.);
                              }

        TH1F *dedx_p2[48];
        for(int i=0;i<48;i++) {
        sprintf(htitle,"dedx_p2_%d",(i+1)); 
        dedx_p2[i] = new TH1F(htitle,htitle,100,0.,5.);
                              }

        TH1F *dedx_he2[48];
        for(int i=0;i<48;i++) {
        sprintf(htitle,"dedx_he2_%d",(i+1)); 
        dedx_he2[i] = new TH1F(htitle,htitle,100,0.,15.);
                              }

        TH1F *dedx_c2[48];
        for(int i=0;i<48;i++) {
        sprintf(htitle,"dedx_c2_%d",(i+1)); 
        dedx_c2[i] = new TH1F(htitle,htitle,100,10.,100.);
                              }


        TH1F *invbeta[164];
        for(int i=0;i<164;i++) {
        sprintf(htitle, "invbeta_%d",(i+1));
        invbeta[i] = new TH1F(htitle,htitle,100,0.2,1.8);
                              }

        TH1F *pos[24];
        for(int i=0;i<24;i++) {
        sprintf(htitle, "pos_%d",(i+1));
        pos[i] = new TH1F(htitle,htitle,100,-10.,10.);
                              }



//------------------------------------------------------------------


   //   Create a timer object to benchmark this loop
    TStopwatch timer;
    timer.Start();

     Float_t fRealTime=0.;
     Double_t cputime,rtime,rate;


   TString wd=gSystem->WorkingDirectory();
    //
    //  Elena  Ischia talk

    PamLevel2*  pam_event  = new PamLevel2(dir,list,"+AUTO");       // << create pamela event

//    CaloNuclei *cn = new CaloNuclei(pam_event);
    CaloNuclei *cn = new CaloNuclei(pam_event,alg);

//    ToFNaNuclei *tna = new ToFNaNuclei();
//    tna->InitPar("/gpfs/wizard/flight/analysis/menn/devel/9th/tof/flight/ToFNaNuclei/GoodParam/");


    ToFdEdx_patch *tdedx = new ToFdEdx_patch();
    tdedx->InitPar("/gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFdEdx_patch/GoodParam/","flight");

    ToFdEdx_patch *tdedx1 = new ToFdEdx_patch();
    tdedx1->InitPar("/gpfs/wizard/flight/analysis/menn/devel/10th/tof/flight/ToFdEdx_patch/GoodParam/","flight");

//    ToFPatch *tofpatch = new ToFPatch();
//    tofpatch->InitPar("/gpfs/wizard/flight/analysis/menn/devel/9th/tof/flight/ToFPatch_v5/GoodParam/");


Int_t ipad_a[6] = {0,8,14,16,18,21};
Int_t ipmt_a[6] = {0,16,28,32,36,42};

   
    Int_t nevent = pam_event->GetEntries();

//--------------------------------------------------------------------

    Int_t nbins = 1500;

    Double_t  thelp;

        for(int i=0;i<48;i++) {
        sprintf(htitle,"sc_adc_%d",(i+1)); 
        sc_adc_he[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,0.,2000.);

        sprintf(htitle,"sc_dedx_p1_%d",(i+1)); 
        sc_dedx_p1[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,0.,5.);

        sprintf(htitle,"sc_dedx_he1_%d",(i+1)); 
        sc_dedx_he1[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,0.,15.);

        sprintf(htitle,"sc_dedx_c1_%d",(i+1)); 
        sc_dedx_c1[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,10.,80.);

        sprintf(htitle,"sc_dedx_p2_%d",(i+1)); 
        sc_dedx_p2[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,0.,5.);

        sprintf(htitle,"sc_dedx_he2_%d",(i+1)); 
        sc_dedx_he2[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,0.,15.);

        sprintf(htitle,"sc_dedx_c2_%d",(i+1)); 
        sc_dedx_c2[i] = new TH2F(htitle,htitle,nbins,tmin,tmax,100,10.,80.);
                              }


       TH2F *trkdedx_rig1 = new TH2F("trkdedx_rig1","",100,0.1,10,100,0.,50);
       TH2F *trkdedx_rig2 = new TH2F("trkdedx_rig2","",100,0.1,10,100,0.,50);
       TH2F *zcalo_rig   = new TH2F("zcalo_rig","",100,0.5,10,100,0.,10);
       TH2F *zcalo_rig1   = new TH2F("zcalo_rig1","",100,0.5,10,100,0.,10);

//--------------------------------------------------------------------

    Int_t totevent = -1;
    Int_t trkevent = -1;
    Int_t p_event = -1;
    Int_t he_event = -1;
    Int_t c_event = -1;
    Int_t invb_event = -1;
 
    Int_t paddlestrack[6]={0};
 
    Int_t trigconfsaved = 0;

    Double_t absolute_time;
    UInt_t atime;

    Double_t  tsum = 0.;
    Double_t  nsum = 0.;

    Double_t  yhelp,yhelp0,yhelp1,yhelp2,yhelp3;
    Double_t adcl,adcr,adclo,adcro,am,theta,xkorr,dx,dy,dr;
    Int_t ihelp,jj,kk,hh;

    Float_t na_dedx[48];  // ToFNaNuclei dEdx
    Float_t patch_dedx[48];  // ToFPatch dEdx
    Float_t dedx_wm[48];  // new ToFdEdx_patch
    Float_t dedx_wm1[48];  // new ToFdEdx_patch


    //

    Int_t ntrk1=0;
    Int_t ntrk2=0;
    if(nevent < nmax)nmax=nevent;
    cout << endl<<" Start loop over "<<nmax<<" events   "<<endl;

    for (Int_t iev=0; iev<nmax;iev++){
        //

        pam_event->Clear();
        pam_event->GetEntry(iev);

//----- Define time
       thelp=pam_event->GetOrbitalInfo()->absTime;
       atime = thelp;

       if (iev==0) tmin = thelp;  // first event

   if ( iev == 0) cout<<"First event time : "<<thelp<<endl;   
//   if ( ntrk1%50000 == 0) cout<<iev<<" "<<thelp<<" ntrk  "<<ntrk1<<" H "<<p_event<<" He "<<he_event<<" C "<<c_event<<" 1/beta "<<invb_event<<endl;   
   if ( iev%50000 == 0) {
        timer.Stop();
        fRealTime = fRealTime + timer.RealTime();
        timer.Start();
        rate=(1.*iev)/fRealTime;
       cout<<" RealTime= "<<fRealTime<<" Rate "<<rate<<" events/sec"<<endl;
       cout<<iev<<" "<<thelp<<" ntrk  "<<ntrk1<<" H "<<p_event<<" He "<<he_event<<" C "<<c_event<<" 1/beta "<<invb_event<<endl;   
                          }



      Int_t ip = 0;
      Int_t ihe = 0;
      Int_t ic=0;
      
      
  Int_t istandard = 1;

  if (istandard == 1) {


	//////////////////////////////////////////////////////////
	///
	/// standard 
	///
	//////////////////////////////////////////////////////////


// old
//    int npt = pam_event->GetTrkLevel2()->GetNTracks();
// new
	  int npt = pam_event->GetNTracks(alg); 

    if( npt>0 ) ntrk1++;  //if n.physical tracks>0

//=======================  only  single tracks    ========================
  
    if( npt>0 ) { //if n.physical tracks>0

/*
    PamTrack *pamtrack = pam_event->GetTrack(0);       //this method selects the 0th physical track, discarding the image	   
    TrkTrack* trktrack = pamtrack->GetTrkTrack();
    ToFTrkVar* toftrack = pamtrack->GetToFTrack(); 
    CaloTrkVar* calotrack = pamtrack->GetCaloTrack(); 
*/

    PamTrack *pamtrack = pam_event->GetTrack(0,alg); //the method selects the 0th physical track, discarding the image	
    ExtTrack* trktrack = pamtrack->GetExtTrack();
    ToFTrkVar* toftrack = pamtrack->GetToFTrack(); 
    CaloTrkVar *calotrack = pam_event->GetTrack(0,alg)->GetCaloTrack();





     if(
       npt ==1 && //the tracks is single
       trktrack->chi2 > 0 && //the fit converged
    true){


Float_t chi=trktrack->chi2;
Float_t deflection=trktrack->al[4];
Float_t rig = 1./deflection;
Float_t beta =  toftrack->beta[12] ;

//cout<<"beta "<<beta<<endl;

//cout<<atime<<" "<<rig<<" "<<"beta "<<beta<<" chi "<<endl;

            if(
                chi > 0    &&
                pow((double)chi,0.25) < (3.5+1.85*fabs(deflection)) &&
                trktrack->GetNX() >= 4  &&
                trktrack->GetNY() >= 3  &&
                deflection > 0 &&
                deflection < 10 &&
//                beta > 0.1 &&
//                beta < 2.0 &&
//                beta > 0.1 &&
//                beta < 0.6 &&
                true
                ){


//cout<<"passing trk cuts "<<endl;
//gets(s);

//cout<<iev<<"  STD rig  "<<1/deflection<<" chi "<<chi<<"  beta "<<beta<<endl;

//-----------------------------------------------------------------------



      Double_t def = trktrack->al[4];
      yhelp0 = 0.5 + 0.6*def*def ;
      yhelp1 = 4. + 0.8*def*def ;
      yhelp2 = 4. + 8.*def*def;
      yhelp3 = 10. + 12.5*def*def;

      yhelp = trktrack->GetDEDX(); 

      trkdedx_rig1->Fill((1./def),(yhelp/1.3));

if ((yhelp>yhelp0) && (yhelp<yhelp1)) ip=1;   //protons
if ((yhelp>yhelp2) && (yhelp<yhelp3)) ihe=1;  //helium
if (yhelp>yhelp3) ic=1;  // Z>2

//cout<<" rig "<<rig<<"  trktrack->GetDEDX() "<<yhelp<<" ip= "<<ip<<endl;

if (ip==1) p_event = p_event + 1 ;
if (ihe==1) he_event = he_event + 1 ;

   if (ihe==1) dedx_he_trk->Fill(yhelp/1.3);


Float_t xv[6],yv[6];
for (jj=0; jj<6; jj++){
xv[jj]=toftrack->xtr_tof[jj];
yv[jj]=toftrack->ytr_tof[jj];
}

Float_t xt,yt,xl,xh,yl,yh;
//Float_t beta;
Int_t npad,i1,i2;

Int_t paddleidoftrack[6]= {0} ;

for (Int_t  ilay=0; ilay<6; ilay++) {
paddleidoftrack[ilay] = pam_event->GetToFLevel2()->GetPaddleIdOfTrack(xv[ilay], yv[ilay], ilay, 0.0) ;
                                    }

Int_t  tof11_i,tof12_i,tof21_i,tof22_i,tof31_i,tof32_i;

tof11_i = paddleidoftrack[0];
tof12_i = paddleidoftrack[1];
tof21_i = paddleidoftrack[2];
tof22_i = paddleidoftrack[3];
tof31_i = paddleidoftrack[4];
tof32_i = paddleidoftrack[5];

//cout<<tof11_i<<" "<<tof12_i<<" "<<tof21_i<<" "<<tof22_i<<" "<<tof31_i<<" "<<tof32_i<<endl;

//----------------------------------------------------------------
//----------------------------------------------------------------

Float_t  xhelp;

       jj = tof11_i ;
       if (jj != -1) {
         yhelp = toftrack->ytofpos[0];
         if (fabs(yhelp-yv[0]) > 0.01) pos[0+jj]->Fill(yhelp-yv[0]);
            }
       jj = tof12_i ;
       if (jj != -1) {
         xhelp = toftrack->xtofpos[0];
           if (fabs(xhelp-xv[1]) > 0.01) pos[8+jj]->Fill(xhelp-xv[1]);
            }

       jj = tof21_i ;
       if (jj != -1) {
         xhelp = toftrack->xtofpos[1];
           if (fabs(xhelp-xv[2]) > 0.01) pos[14+jj]->Fill(xhelp-xv[2]);
            }
       jj = tof22_i ;
       if (jj != -1) {
         yhelp = toftrack->ytofpos[1];
           if (fabs(yhelp-yv[3]) > 0.01) pos[16+jj]->Fill(yhelp-yv[3]);
            }

       jj = tof31_i ;
       if (jj != -1) {
         yhelp = toftrack->ytofpos[2];
           if (fabs(yhelp-yv[4]) > 0.01) pos[18+jj]->Fill(yhelp-yv[4]);
            }
       jj = tof32_i ;
       if (jj != -1) {
         xhelp = toftrack->xtofpos[2];
           if (fabs(xhelp-xv[5]) > 0.01) pos[21+jj]->Fill(xhelp-xv[5]);
            }

//======================================================================
//========  select high rigidity protons to check inverse beta  ========
//======================================================================

if ((ip == 1) && (trktrack->al[4] > 0.) && (trktrack->al[4] < 0.08)) {

invb_event++;

//----------------------------------------------------------------
//---------  S1 - S3 ---------------------------------------------
//----------------------------------------------------------------

//---------  S11 - S31 -------------------------------------------

       i1 = tof11_i;
       i2 = tof31_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=i1*3+i2;
         beta =  toftrack->beta[0] ;
         invbeta[ihelp]->Fill(1./beta);
           }

//---------  S11 - S32 -------------------------------------------

       i1 = tof11_i;
       i2 = tof32_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=24+i1*3+i2;
         beta =  toftrack->beta[1] ;
         invbeta[ihelp]->Fill(1./beta);
           }


//---------  S12 - S31 -------------------------------------------

       i1 = tof12_i;
       i2 = tof31_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=48+i1*3+i2;
         beta =  toftrack->beta[2] ;
         invbeta[ihelp]->Fill(1./beta);
           }

//---------  S12 - S32 -------------------------------------------

       i1 = tof12_i;
       i2 = tof32_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=66+i1*3+i2;
         beta =  toftrack->beta[3] ;
         invbeta[ihelp]->Fill(1./beta);
           }


//----------------------------------------------------------------
//---------  S2 - S3 ---------------------------------------------
//----------------------------------------------------------------

//---------  S21 - S31 -------------------------------------------

       i1 = tof21_i;
       i2 = tof31_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=84+i1*3+i2;
         beta =  toftrack->beta[4] ;
         invbeta[ihelp]->Fill(1./beta);
             }

//---------  S21 - S32 -------------------------------------------

       i1 = tof21_i;
       i2 = tof32_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=90+i1*3+i2;
         beta =  toftrack->beta[5] ;
         invbeta[ihelp]->Fill(1./beta);
           }


//---------  S22 - S31 -------------------------------------------

       i1 = tof22_i;
       i2 = tof31_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=96+i1*3+i2;
         beta =  toftrack->beta[6] ;
         invbeta[ihelp]->Fill(1./beta);
               }

//---------  S22 - S32 -------------------------------------------

       i1 = tof22_i;
       i2 = tof32_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=102+i1*3+i2;
         beta =  toftrack->beta[7] ;
         invbeta[ihelp]->Fill(1./beta);
           }


//----------------------------------------------------------------
//---------  S1 - S2 ---------------------------------------------
//----------------------------------------------------------------

//---------  S11 - S21 -------------------------------------------

       i1 = tof11_i;
       i2 = tof21_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=108+i1*2+i2;
         beta =  toftrack->beta[8] ;
         invbeta[ihelp]->Fill(1./beta);
           }


//---------  S11 - S22 -------------------------------------------

       i1 = tof11_i;
       i2 = tof22_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=124+i1*2+i2;
         beta =  toftrack->beta[9] ;
         invbeta[ihelp]->Fill(1./beta);
           }


//---------  S12 - S21 -------------------------------------------

       i1 = tof12_i;
       i2 = tof21_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=140+i1*2+i2;
         beta =  toftrack->beta[10] ;
         invbeta[ihelp]->Fill(1./beta);
           }


//---------  S12 - S22 -------------------------------------------

       i1 = tof12_i;
       i2 = tof22_i;
       if ((i1 != -1) && (i2 != -1)) {
         ihelp=152+i1*2+i2;
         beta =  toftrack->beta[11] ;
         invbeta[ihelp]->Fill(1./beta);
           }

//-------------------------------------------------------------------------------

        }; // charge and R>12


//=======================================================================
//=====================      ToFNaNuclei  ===============================
//=======================================================================
/*
   tna->Process(pam_event,0);
   for (Int_t ipmt=0; ipmt<48; ipmt++) na_dedx[ipmt] = tna->GetdEdx_pmt(ipmt);
*/
//=======================================================================
//=====================      ToFdEdx_patch     ===============================
//=======================================================================

//   tdedx->Process(pam_event,0);
   tdedx->Process(pam_event,alg,tri_or_bi);
   for (Int_t ipmt=0; ipmt<48; ipmt++) dedx_wm[ipmt] = tdedx->GetdEdx_pmt(ipmt);

//=======================================================================
//=====================      ToFPatch     ===============================
//=======================================================================
//   tofpatch->Process(pam_event,0);
//   for (Int_t ipmt=0; ipmt<48; ipmt++) patch_dedx[ipmt] = tofpatch->GetdEdx_pmt(ipmt);
//======================================================================
//===========================  select helium ===========================
//======================================================================

if ((ihe==1) && (deflection<0.5)) {
//if ((ihe==1)) {


// fill adc "raw" values
// we have raw ADC values in then "PMT" class (track independent)

     Float_t adc[48]; //,dedx[48];
     for (jj=0;jj<48;jj++)   adc[jj]=4095.;
//     for (jj=0;jj<48;jj++)   dedx[jj]=1000.;
                

Float_t aa;
Int_t ipm,ich,ihb;


for (Int_t ipmt=0; ipmt<pam_event->GetToFLevel2()->npmt() ; ipmt++){
ToFPMT *tofpmt = pam_event->GetToFLevel2()->GetToFPMT(ipmt);
ipm = tofpmt->pmt_id;
aa = (tofpmt->adc) ;
adc_he[ipm]->Fill(aa);
sc_adc_he[ipm]->Fill(thelp,aa);
adc[ipm] = aa;  // raw ADC values
//cout<<ipm<<" "<<adc[ipm]<<endl;
                 }


//   if (dedx_wm[9]<1000) cout<<adc[9]<<"  "<<dedx_wm[9]<<endl;

//  normal loop for 10th red data

 for (Int_t ipmt=0; ipmt<toftrack->npmtadc; ipmt++){
//          Float_t dEdx = toftrack->dedx[ipmt];
          Int_t pmtadc = toftrack->pmtadc[ipmt];
          Int_t adcflag = toftrack->adcflag[ipmt];

   Int_t iverbose=0;
//   if (pmtadc==14) iverbose=1;
//   if (pmtadc==27) iverbose=1;
//   if (pmtadc==18) iverbose=1;

          Float_t dEdx1 = toftrack->dedx[ipmt];
//          Float_t dEdx1 = na_dedx[pmtadc];  // ToFNaNuclei dEdx
          if (adcflag==0) dedx_he1[pmtadc]->Fill(dEdx1);
          if (adcflag==0) sc_dedx_he1[pmtadc]->Fill(thelp,dEdx1);

                                                     }



//  new loop over hitted paddles for tofdedx_patch data


for (Int_t ilay=0; ilay<6; ilay++) {
Int_t i1 = paddleidoftrack[ilay];
if (i1 != -1) {
 Int_t pmtadc = ipmt_a[ilay]+(i1*2); 

   Int_t iverbose=0;
//   if (pmtadc==14) iverbose=1;
//   if (pmtadc==27) iverbose=1;
//   if (pmtadc==17) iverbose=1;

    Float_t dEdx2 = dedx_wm[pmtadc];  // ToFdEdx_patch  wm

   if (iverbose==1) cout<<atime<<"  "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
    if (adc[pmtadc] < 4095)   dedx_he2[pmtadc]->Fill(dEdx2);
    if (adc[pmtadc] < 4095)   sc_dedx_he2[pmtadc]->Fill(thelp,dEdx2);

   pmtadc = pmtadc + 1;
    dEdx2 = dedx_wm[pmtadc];  // ToFdEdx_patch  wm

    iverbose=0;
//   if (pmtadc==14) iverbose=1;
//   if (pmtadc==27) iverbose=1;
//   if (pmtadc==17) iverbose=1;

   if (iverbose==1) cout<<atime<<"  "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
    if (adc[pmtadc] < 4095)   dedx_he2[pmtadc]->Fill(dEdx2);
    if (adc[pmtadc] < 4095)   sc_dedx_he2[pmtadc]->Fill(thelp,dEdx2);


} // i1>-1
} // ilay



// End fill adc 

  }; //ihe==1


//======================================================================
//======================================================================
//===========   select Protons with Tracker  ====================
//======================================================================

//if ( (deflection<0.5) && (ip == 1) && (fabs(charge-1.)<0.15) ) { // iz = 1 
if ( (deflection<0.5) && (ip == 1) ) { // iz = 1 

//cout<<"Proton detected !!!! "<<endl;


Float_t adc[48];
for (jj=0;jj<48;jj++)   adc[jj]=4095.;

Float_t aa;
Int_t ipm,ich,ihb;


for (Int_t ipmt=0; ipmt<pam_event->GetToFLevel2()->npmt() ; ipmt++){
ToFPMT *tofpmt = pam_event->GetToFLevel2()->GetToFPMT(ipmt);
ipm = tofpmt->pmt_id;
aa = (tofpmt->adc) ;
adc[ipm] = aa;  // raw ADC values
                 }


/*
 for (Int_t ipmt=0; ipmt<toftrack->npmtadc; ipmt++){
  Int_t ii = toftrack->pmtadc[ipmt];
  cout<<"In monitor400  "<<ipmt<<" "<<ii<<" "<<adc[ii]<<endl;                                 
                   }
*/

// normal loop for 10th reduction data

 for (Int_t ipmt=0; ipmt<toftrack->npmtadc; ipmt++){
          Int_t pmtadc = toftrack->pmtadc[ipmt];

   Int_t iverbose=0;
//   if (pmtadc==27) iverbose=1;
//   if (pmtadc==17) iverbose=1;

          Int_t adcflag = toftrack->adcflag[ipmt];

          Float_t dEdx1 = toftrack->dedx[ipmt];
//          Float_t dEdx1 = na_dedx[pmtadc];  // ToFNaNuclei dEdx

         if (iverbose==1) cout<<"10th "<<atime<<"  "<<rig<<" adc "<<adc[pmtadc]<<"  10th dedx "<<dEdx1<<endl;
   
          if (adcflag==0) dedx_p1[pmtadc]->Fill(dEdx1);
          if (adcflag==0) sc_dedx_p1[pmtadc]->Fill(thelp,dEdx1);

                                                      
if (iverbose==1) gets(s);
                                              }  // ipmt
                                              
                                              
 //  new loop over hitted paddles for tofdedx_patch data


for (Int_t ilay=0; ilay<6; ilay++) {
Int_t i1 = paddleidoftrack[ilay];
if (i1 != -1) {
 Int_t pmtadc = ipmt_a[ilay]+(i1*2); 

   Int_t iverbose=0;
//   if (pmtadc==14) iverbose=1;
//   if (pmtadc==27) iverbose=1;
//   if (pmtadc==17) iverbose=1;

    Float_t dEdx2 = dedx_wm[pmtadc];  // ToFdEdx_patch  wm

   if (iverbose==1) cout<<"tofdedx_patch "<<atime<<"  "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
    if (adc[pmtadc] < 4095)   dedx_p2[pmtadc]->Fill(dEdx2);
    if (adc[pmtadc] < 4095)   sc_dedx_p2[pmtadc]->Fill(thelp,dEdx2);

if (iverbose==1) gets(s);

   pmtadc = pmtadc + 1;
    dEdx2 = dedx_wm[pmtadc];  // ToFdEdx_patch  wm

    iverbose=0;
//   if (pmtadc==14) iverbose=1;
//   if (pmtadc==27) iverbose=1;
//  if (pmtadc==17) iverbose=1;

   if (iverbose==1) cout<<"tofdedx_patch "<<atime<<"  "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
    if (adc[pmtadc] < 4095)   dedx_p2[pmtadc]->Fill(dEdx2);
    if (adc[pmtadc] < 4095)   sc_dedx_p2[pmtadc]->Fill(thelp,dEdx2);

if (iverbose==1) gets(s);

} // i1>-1
} // ilay

                                             
                                              
                                              
               }  // protons selected 

//======================================================================

          };  // tracking cut
       };  // fit converged
      };   //  ntrk = 1 or 2
      
//======================================================================


} // istandard



Int_t inuc = 1;

if (inuc == 1) {

	//////////////////////////////////////////////////////////
	///
	/// NUC 
	///
	//////////////////////////////////////////////////////////

	  int npt = pam_event->GetNTracks(alg1); 

    if( npt>0 ) ntrk2++;  //if n.physical tracks>0

//=======================  only  single tracks    ========================
  
    if( npt>0 ) { //if n.physical tracks>0

    PamTrack *pamtrack1 = pam_event->GetTrack(0,alg1); //the method selects the 0th physical track, discarding the image	

    ExtTrack* trktrack1 = pamtrack1->GetExtTrack();
    ToFTrkVar* toftrack1 = pamtrack1->GetToFTrack(); 
    CaloTrkVar *calotrack1 = pam_event->GetTrack(0,alg1)->GetCaloTrack();


     if(
       npt ==1 && //the tracks is single
        trktrack1->chi2 > 0 && //the fit converged
    true){


            //

Float_t chi=trktrack1->chi2;
Float_t deflection=trktrack1->al[4];
Float_t rig = 1./deflection;
Float_t beta =  toftrack1->beta[12] ;

            if(
                chi > 0    &&
//                pow((double)chi,0.25) < (3.5+1.85*fabs(deflection)) &&
                trktrack1->GetNX() >= 3  &&
                trktrack1->GetNY() >= 2  &&
                trktrack1->GetLeverArmX() >= 4 &&
                deflection > 0 &&
                deflection < 10 &&
//                beta > 0.1 &&
//                beta < 0.6 &&
                true
                ){

//cout<<iev<<"  NUC rig  "<<1/deflection<<" chi "<<chi<<"  beta "<<beta<<endl;


      Double_t def = trktrack1->al[4];
      yhelp0 = 0.5 + 0.6*def*def ;
      yhelp1 = 4. + 0.8*def*def ;
      yhelp2 = 4. + 8.*def*def;
      yhelp3 = 10. + 12.5*def*def;

      yhelp = trktrack1->GetDEDX(); 

//      trkdedx_rig1->Fill((1./def),(yhelp/1.3));

if ((yhelp>yhelp0) && (yhelp<yhelp1)) ip=1;   //protons
if ((yhelp>yhelp2) && (yhelp<yhelp3)) ihe=1;  //helium
if (yhelp>yhelp3) ic=1;  // Z>2

if (ip==1) p_event = p_event + 1 ;
if (ihe==1) he_event = he_event + 1 ;

//   if (ihe==1) dedx_he_trk->Fill(yhelp/1.3);




Float_t xv[6],yv[6];
for (jj=0; jj<6; jj++){
xv[jj]=toftrack1->xtr_tof[jj];
yv[jj]=toftrack1->ytr_tof[jj];
}

Int_t paddleidoftrack[6]= {0} ;

for (Int_t  ilay=0; ilay<6; ilay++) {
paddleidoftrack[ilay] = pam_event->GetToFLevel2()->GetPaddleIdOfTrack(xv[ilay], yv[ilay], ilay, 0.0) ;
                                    }


//======================================================================
//=================   select Carbon with Calo   ========================
//======================================================================


Float_t charge = 1000.;
Float_t mip=0;

//Float_t charge1 = 1000;
//charge1  = cn->Get_charge_siegen1();

// mip=cn->Get_StdEdx1();


charge=1000;

// Data from file Calo_Bands_New_7.dat
Float_t C0[9] = {0 , 1 , 2 , 3 , 4 , 5 , 6 , 8 , 90 };
Float_t B0[9] = {0 , -2.03769 , 7.61781 , 19.7098 , 60.5598 , 57.9226 , 14.8368 , -1358.83 , 8200 };
Float_t B1[9] = {0 , 0.0211274 , 9.32057e-010 , 4.47241e-07 , 1.44826e-06 , 2.6189e-05 , 0.00278178 , 55.5445 , 0 };
Float_t B2[9] = {0 , -3.91742 , -20.0359 , -16.3043 , -16.9471 , -14.4622 , -10.9594 , -2.38014 , 0 };
Float_t B3[9] = {0 , 11.1469 , -6.63105 , -27.8834 , -132.044 , -55.341 , 173.25 , 4115 , 0 };
Float_t B4[9] = {0 , -14.3465 , -0.485215 , 18.8122 , 117.533 , -14.0898 , -325.269 , -4388.89 , 0 };
Float_t B5[9] = {0 , 6.24281 , 3.96018 , 0 , -26.1881 , 42.9731 , 182.697 , 1661.01 , 0 };


Float_t x1[9],y1[9];
Int_t n1 = 9;

//    if ( !calotrack ) printf(" ERROR: No CaloTrack!\n");

  Int_t view = 0;
  Int_t plane = 0;
  Int_t strip = 0;

  Int_t indx = 0;
  Float_t vfpl[96];
  Int_t stfpl[96];
  memset(vfpl, 0, 96*sizeof(Float_t));
  memset(stfpl, 0, 96*sizeof(Int_t));

  Float_t  dedx1 = 0.;
  Float_t  stdedx1 = 0.;
  Float_t  maxrel = 0.;
  
  mip=0.;
  
  for ( Int_t is=0; is<pam_event->GetCaloLevel1()->istrip; is++ ){
    //
    mip = pam_event->GetCaloLevel1()->DecodeEstrip(is,view,plane,strip);

      if ( strip != -1 && 
	   view == 1   && 
	   plane == 0  &&
	   ( strip == (calotrack1->tibar[0][1]-1) || strip == (calotrack1->tibar[0][1]-2) || strip == (calotrack1->tibar[0][1]) ) 
	   && true ){      
	   dedx1 += mip;
      };
      
    if ( strip != -1 && view == 1 && plane == 0 ) {
      stfpl[indx] = strip;
      vfpl[indx] = mip;
      indx++;
    };

      
}  // is


  // find energy released along the strip of maximum on the first plane and on the two neighbour strips
  //
  if ( indx > 0 ){
    Int_t mindx = (Int_t)TMath::LocMax(indx,vfpl);
    for (Int_t ii=0; ii<indx; ii++){
      if ( stfpl[ii] == stfpl[mindx] ) stdedx1 += vfpl[ii];
      if ( (mindx-1)>=0 && stfpl[ii] == (stfpl[mindx]-1) ) stdedx1 += vfpl[ii];
      if ( (mindx+1)<96 && stfpl[ii] == (stfpl[mindx]+1) ) stdedx1 += vfpl[ii];
      //      if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii];
      //      if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii];
    };
  maxrel = vfpl[mindx];
  } else {
    stdedx1 = 0.;
    maxrel = 0.;
  };
  
  
  // cout<<"cn->Get_stdedx1 "<<mip1<<" selfmade dedx1 "<<dedx1<<" selfmade stdedx1 "<<stdedx1<<"\n";
  //
  //cout<<"Calo dedx1 "<<dedx1<<endl;

//-----------------------------------------------------------------------


    beta =  toftrack1->beta[12] ;

//    if (beta<2.) {  // it makes no sense to allow beta=5 or so...
    if ((beta>0.1)&&(beta<2.)) {  // it makes no sense to allow beta=5 or so...  reject negative beta

//------ Calo function

    Float_t  mip;
//    mip=dedx1;
//   mip=cn->Get_StdEdx1();
     mip = stdedx1 ;
    

    if (mip>0)   {

    Float_t betahelp =   pow(beta, 1.8);
    Float_t ym = mip*betahelp;
    Float_t xb = beta;

    for ( Int_t jj=0; jj<9; jj++ ){
    x1[jj] = B0[jj]+B1[jj]*pow(xb,B2[jj])+B3[jj]*xb+B4[jj]*xb*xb+B5[jj]*xb*xb*xb;
    y1[jj] = C0[jj]*C0[jj] ;
                                  }

    TGraph *gr1 = new TGraph(n1,x1,y1);
    TSpline3 *spl1 = new TSpline3("grs",gr1);    // use a cubic spline
    Float_t chelp = spl1->Eval(ym);
    charge = TMath::Sqrt(chelp);
    gr1->Delete();
    spl1->Delete();

                    } // if (mip1>0)
                    } // beta < 2

//cout<<"Handmade charge "<<charge<<endl;
//cout<<"Handmade : stdedx1 "<<stdedx1<<" charge "<<charge<<endl;

//cout<<"Handmade charge "<<charge<<" new CaloNuclei "<<charge1<<endl;


//cout<<iev<<" rig "<<(1./def)<<" beta "<<beta<<" stdedx1 "<<stdedx1<<" charge "<<charge<<endl;

//=================================================================

//cout<<"CaloNuclei : stdedx1 "<<mip<<" charge "<<charge<<endl;

zcalo_rig->Fill((1./def),charge);

if (ic==1) { // otherwise  the count is swamped by background of  protons...

 zcalo_rig1->Fill((1./def),charge);

ic=0;
//if (fabs(charge-6.)<0.40) ic = 1;
if ((5.6<charge)&&(charge<6.4)) ic = 1;

if (ic==1) c_event = c_event + 1 ;

//=======================================================================
//=======================================================================
//=======================================================================

if (fabs(charge-6.)<0.40) { // iz = 6 
if (deflection<0.5) {


//	cout<<" main program:  npt = pam_event->GetNTracks(alg1) "<<npt<<endl;

//=====================      ToFdEdx_patch     ==========================

   tdedx1->Process(pam_event,alg1);
   for (Int_t ipmt=0; ipmt<48; ipmt++) dedx_wm1[ipmt] = tdedx1->GetdEdx_pmt(ipmt);

//=======================================================================

Float_t adc[48];
for (jj=0;jj<48;jj++)   adc[jj]=4095.;

Float_t aa;
Int_t ipm,ich,ihb;


for (Int_t ipmt=0; ipmt<pam_event->GetToFLevel2()->npmt() ; ipmt++){
ToFPMT *tofpmt = pam_event->GetToFLevel2()->GetToFPMT(ipmt);
ipm = tofpmt->pmt_id;
aa = (tofpmt->adc) ;
adc[ipm] = aa;  // raw ADC values
                 }



 for (Int_t ipmt=0; ipmt<toftrack1->npmtadc; ipmt++){
          Int_t pmtadc = toftrack1->pmtadc[ipmt];
//          Float_t dEdx = toftrack1->dedx[ipmt];
//          Float_t dEdx10th = toftrack->dedx[ipmt];
//          cout<<ipmt<<" 10th "<<dEdx10th<<" TNa "<<dEdx<<endl;
          Int_t adcflag = toftrack1->adcflag[ipmt];

          Float_t dEdx1 = toftrack1->dedx[ipmt];
//         Float_t dEdx1 = na_dedx[pmtadc];  // ToFNaNuclei dEdx

   Int_t iverbose=0;
//   if (pmtadc==33) iverbose=1;
   if (iverbose==1) cout<<ipmt<<" "<<pmtadc<<" adc "<<adc[pmtadc]<<" 10th "<<dEdx1<<endl;
 
          if (adcflag==0) dedx_c1[pmtadc]->Fill(dEdx1);
          if (adcflag==0) sc_dedx_c1[pmtadc]->Fill(thelp,dEdx1);

                                                     }


                                              
 //  new loop over hitted paddles for tofdedx_patch data


for (Int_t ilay=0; ilay<6; ilay++) {
Int_t i1 = paddleidoftrack[ilay];
if (i1 != -1) {
 Int_t pmtadc = ipmt_a[ilay]+(i1*2); 

   Int_t iverbose=0;
//   if (pmtadc==14) iverbose=1;
//   if (pmtadc==27) iverbose=1;
//   if (pmtadc==33) iverbose=1;

    Float_t dEdx2 = dedx_wm1[pmtadc];  // ToFdEdx_patch  wm

//   if (iverbose==1) cout<<atime<<"  "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
  if (iverbose==1) cout<<pmtadc<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
    if (adc[pmtadc] < 4095)   dedx_c2[pmtadc]->Fill(dEdx2);
    if (adc[pmtadc] < 4095)   sc_dedx_c2[pmtadc]->Fill(thelp,dEdx2);

   pmtadc = pmtadc + 1;
    dEdx2 = dedx_wm1[pmtadc];  // ToFdEdx_patch  wm

    iverbose=0;
//   if (pmtadc==14) iverbose=1;
//   if (pmtadc==27) iverbose=1;
//   if (pmtadc==33) iverbose=1;

 //   cout<<pmtadc<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
  if (iverbose==1) cout<<pmtadc<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;

 //  if (iverbose==1) cout<<atime<<"  "<<rig<<" adc "<<adc[pmtadc]<<" tofdedxpatch "<<dEdx2<<endl;
    if (adc[pmtadc] < 4095)   dedx_c2[pmtadc]->Fill(dEdx2);
    if (adc[pmtadc] < 4095)   sc_dedx_c2[pmtadc]->Fill(thelp,dEdx2);


} // i1>-1
} // ilay


               } // deflection < 0.5
               } //  charge-6  < 0.4
               
//=======================================================================
//=======================================================================
//=======================================================================

   }



//======================================================================

          };  // tracking cut
       };  // fit converged
      };   //  ntrk = 1 or 2
      
//======================================================================

} // inuc

//======================================================================




thelp=pam_event->GetOrbitalInfo()->absTime;

tsum = tsum + thelp ;
nsum = nsum + 1. ;


//---------------------------------------------------------------------      
//---------------------------------------------------------------------        
     
    };    // end loop over the events

//---------------------------------------------------------------------      
//---------------------------------------------------------------------      

    cout << endl << endl << " Done "<< endl<<endl;
    cout <<"Std " << ntrk1 <<" tracks over "<<nmax<<" events ("<< 100*ntrk1/nmax<<"%)"<<endl;
    cout <<"NUC " << ntrk2 <<" tracks over "<<nmax<<" events ("<< 100*ntrk2/nmax<<"%)"<<endl;
    cout<<" H "<<p_event<<" He "<<he_event<<" C "<<c_event<<" 1/beta "<<invb_event<<endl;   


//  Stop timer and print results
    timer.Stop();
    rtime = timer.RealTime();
    cputime = timer.CpuTime();
    printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,cputime);


//================================================================================
//================================================================================

Double_t xhelp = tsum/nsum ;  // mean time in interval


cout<<" mean time in interval "<<xhelp<<endl;



 Float_t XX[1],YY[1],eXX[1],eYY[1];

     XX[0]  = xhelp;
     eXX[0] = 0.00001;

// TGraphErrors *gr_ibeta[164];
// TGraphErrors *gr_pos[24];

 TGraphErrors *gr_adc[48];
 
 TGraphErrors *gr_dedx_p1[48];
 TGraphErrors *gr_dedx_he1[48];
 TGraphErrors *gr_dedx_c1[48];

 TGraphErrors *gr_dedx_p2[48];
 TGraphErrors *gr_dedx_he2[48];
 TGraphErrors *gr_dedx_c2[48];
 


TF1 *f1 = new TF1("f1","landau",0,100);
//TF1 *f1 = new TF1("f1","gaus",0.,100.);

   for (jj=0;jj<48;jj++) {

     Float_t xhelp1  = adc_he[jj]->GetMean();
     Float_t xhelp2 = 0.01;
     YY[0]  = xhelp1;
     eYY[0] = xhelp2;
     gr_adc[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);

     f1->SetParameters(0.2*dedx_p1[jj]->GetEntries(), dedx_p1[jj]->GetMean(), 0.5*dedx_p1[jj]->GetRMS());
     f1->SetRange(0.5*dedx_p1[jj]->GetMean(), 2.0*dedx_p1[jj]->GetMean());
     dedx_p1[jj]->Fit("f1", "R");
     Float_t xhelp3 = f1->GetParameter(1);
     Float_t xhelp4 = f1->GetParError(1);
     YY[0]  = xhelp3;
     eYY[0] = xhelp4;
     gr_dedx_p1[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);

     f1->SetParameters(0.2*dedx_he1[jj]->GetEntries(), dedx_he1[jj]->GetMean(), 0.5*dedx_he1[jj]->GetRMS());
     f1->SetRange(0.5*dedx_he1[jj]->GetMean(), 2.0*dedx_he1[jj]->GetMean());
     dedx_he1[jj]->Fit("f1", "R");
     Float_t xhelp5 = f1->GetParameter(1);
     Float_t xhelp6 = f1->GetParError(1);
     YY[0]  = xhelp5;
     eYY[0] = xhelp6;
     gr_dedx_he1[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);

     f1->SetParameters(0.2*dedx_c1[jj]->GetEntries(), dedx_c1[jj]->GetMean(), 0.5*dedx_c1[jj]->GetRMS());
     f1->SetRange(0.5*dedx_c1[jj]->GetMean(), 2.0*dedx_c1[jj]->GetMean());
     dedx_c1[jj]->Fit("f1", "R");
     Float_t xhelp7 = f1->GetParameter(1);
     Float_t xhelp8 = f1->GetParError(1);
     YY[0]  = xhelp7;
     eYY[0] = xhelp8;
     gr_dedx_c1[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);


     f1->SetParameters(0.2*dedx_p2[jj]->GetEntries(), dedx_p2[jj]->GetMean(), 0.5*dedx_p2[jj]->GetRMS());
     f1->SetRange(0.5*dedx_p2[jj]->GetMean(), 2.0*dedx_p2[jj]->GetMean());
     dedx_p2[jj]->Fit("f1", "R");
      xhelp3 = f1->GetParameter(1);
      xhelp4 = f1->GetParError(1);
     YY[0]  = xhelp3;
     eYY[0] = xhelp4;
     gr_dedx_p2[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);

     f1->SetParameters(0.2*dedx_he2[jj]->GetEntries(), dedx_he2[jj]->GetMean(), 0.5*dedx_he2[jj]->GetRMS());
     f1->SetRange(0.5*dedx_he2[jj]->GetMean(), 2.0*dedx_he2[jj]->GetMean());
     dedx_he2[jj]->Fit("f1", "R");
      xhelp5 = f1->GetParameter(1);
      xhelp6 = f1->GetParError(1);
     YY[0]  = xhelp5;
     eYY[0] = xhelp6;
     gr_dedx_he2[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);

     f1->SetParameters(0.2*dedx_c2[jj]->GetEntries(), dedx_c2[jj]->GetMean(), 0.5*dedx_c2[jj]->GetRMS());
     f1->SetRange(0.5*dedx_c2[jj]->GetMean(), 2.0*dedx_c2[jj]->GetMean());
     dedx_c2[jj]->Fit("f1", "R");
      xhelp7 = f1->GetParameter(1);
      xhelp8 = f1->GetParError(1);
     YY[0]  = xhelp7;
     eYY[0] = xhelp8;
     gr_dedx_c2[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);


//  cout<<jj<<" "<<xhelp1<<" "<<xhelp2<<" "<<xhelp3<<" "<<xhelp4<<" "<<xhelp5<<" "<<xhelp6<<" "<<xhelp7<<" "<<xhelp8<<"\n";
// fout1<<jj<<" "<<xhelp1<<" "<<xhelp2<<" "<<xhelp3<<" "<<xhelp4<<" "<<xhelp5<<" "<<xhelp6<<" "<<xhelp7<<" "<<xhelp8<<"\n";
                     }
                     
/*                     
   TF1 *f2 = new TF1("f2","gaus",0.3,1.6);

                                    
     for (jj=0;jj<164;jj++) {
     invbeta[jj]->Fit("f2");
     Float_t xhelp5 = f2->GetParameter(1);
     Float_t xhelp6 = f2->GetParError(1);
     Float_t xhelp7 = f2->GetParameter(2);
     cout<<jj<<" "<<xhelp5<<" "<<xhelp6<<" "<<xhelp7<<"\n";
     YY[0]  = xhelp5;
     eYY[0] = xhelp6;
     gr_ibeta[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);            
                              }

     for (jj=0;jj<24;jj++) {
     pos[jj]->Fit("f2");
     Float_t xhelp8  = f2->GetParameter(1);
     Float_t xhelp9  = f2->GetParError(1);
     Float_t xhelp10 = f2->GetParameter(2);
     cout<<jj<<" "<<xhelp8<<" "<<xhelp9<<" "<<xhelp10<<"\n";
     YY[0]  = xhelp8;
     eYY[0] = xhelp9;
     gr_pos[jj] = new TGraphErrors(1,XX,YY,eXX,eYY);
                         }
*/

//================================================================================
//================================================================================

cout<<"Finished gr-calculation"<<endl;
//cout<<" mean time in interval "<<xhelp<<endl;

  //  -------------------
  //  create output files
  //  -------------------

  TFile *outfh = new TFile(outputname,"RECREATE");

        for(int i=0;i<48;i++) {
        sc_adc_he[i]->Write();
        sc_dedx_p1[i]->Write();
        sc_dedx_he1[i]->Write();
        sc_dedx_c1[i]->Write();
        sc_dedx_p2[i]->Write();
        sc_dedx_he2[i]->Write();
        sc_dedx_c2[i]->Write();
                              }

/*
        for(int i=0;i<164;i++) {
        sprintf(info,"gr_ibeta_%d",i);         
        gr_ibeta[i]->Write(info);
                              }

        for(int i=0;i<24;i++) {
        sprintf(info,"gr_pos_%d",i);         
        gr_pos[i]->Write(info);
                              }
*/

        for(int i=0;i<48;i++) {
        sprintf(info,"gr_adc_%d",i+1);   
 //       cout<<info<<endl;      
        gr_adc[i]->Write(info);
  
        sprintf(info,"gr_dedx_p1_%d",i+1);         
 //       cout<<info<<endl;      
        gr_dedx_p1[i]->Write(info);
        sprintf(info,"gr_dedx_he1_%d",i+1);         
 //       cout<<info<<endl;      
        gr_dedx_he1[i]->Write(info);
        sprintf(info,"gr_dedx_c1_%d",i+1);         
  //      cout<<info<<endl;      
        gr_dedx_c1[i]->Write(info);
        sprintf(info,"gr_dedx_p2_%d",i+1);         
  //      cout<<info<<endl;      
        gr_dedx_p2[i]->Write(info);
        sprintf(info,"gr_dedx_he2_%d",i+1);         
  //      cout<<info<<endl;      
        gr_dedx_he2[i]->Write(info);
        sprintf(info,"gr_dedx_c2_%d",i+1);         
  //      cout<<info<<endl;      
        gr_dedx_c2[i]->Write(info);
 
                              }


       trkdedx_rig1->Write();
       trkdedx_rig2->Write();
       zcalo_rig->Write();
       zcalo_rig1->Write();

  outfh->Close();


// end WM



    gSystem->ChangeDirectory(wd);

};


