#include <CaloNuclei.h>
#include <TGraph.h>
#include <TSpline.h>
#include <TMVA/TSpline1.h>

//--------------------------------------
/**
 * Default constructor 
 */
CaloNuclei::CaloNuclei(){
  Clear();
};

CaloNuclei::CaloNuclei(PamLevel2 *l2p){  
  //
  Clear();
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  N = 5;
  R = 3;
  //
  debug = false;
  // debug = true;
  usetrack = true;
  //
};

void CaloNuclei::Clear(){
  //
  UN = 0;
  tr = 0;
  sntr = 0;
  interplane = 0;
  preq = 0.;
  postq = 0.;
  stdedx1 = 0.;
  ethr = 0.;
  dedx1 = 0.;
  dedx3 = 0.;
  qpremean = 0.; 
  qpremeanN = 0.; 
  maxrel = 0;
  qNmin1 = 0;
  qNmin1_w = 0;
  charge_siegen1 = 0;
  ZCalo_maxrel_b = 0;
  ZCalo_dedx_b = 0;
  ZCalo_dedx_defl= 0;
  ZCalo_Nmin1_defl= 0;
  //
  multhit = false;
  gap = false;
  //
};

void CaloNuclei::Print(){
  //
  Process();
  //
  printf("========================================================================\n");
  printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);
  printf(" interplane [number of available dE/dx before interaction]:....... %i\n",interplane);
  printf(" ethr [threshold used to determine interplane]:................... %f \n",ethr); 
  printf(" dedx1 [dE/dx from the first calorimeter plane]:.................. %f \n",dedx1);
  printf(" stdedx1 [dE/dx from the first calorimeter plane standalone]:..... %f \n",stdedx1);
  printf(" dedx3 [dE/dx (average) if the first 3 Si planes]:................ %f \n",dedx3);
  printf(" multhit [true if interplane determined by multiple hits]:........ %i \n",multhit);
  printf(" gap [true if interplane determined by a gap]:.................... %i \n",gap);
  printf(" preq [total energy in MIP before the interaction plane]:......... %f \n",preq);
  printf(" postq [total energy in MIP after the interaction plane]:......... %f \n",postq);
  printf(" qpremean [truncated mean using 3 planes and 3 strips]:........... %f \n",qpremean); 
  printf(" N [no of used plane]:............................................ %i \n",N); 
  printf(" R [no strip used per plane ]:.................................... %i \n",R); 
  printf(" qpremeanN [truncated mean using N planes and R strips]:.......... %f \n",qpremeanN); 
  printf(" qNmin1 [truncated mean using N-1 planes and R strips]: .......... %f \n",qNmin1); 
  printf(" maxrel [dE/dx of strip with maximum release (I plane)]:.......... %f \n",maxrel); 
  printf(" ZCalo_maxrel_b [Z from maximum release in I Calo plane vs beta].. %f \n",ZCalo_maxrel_b);
  printf(" ZCalo_dedx_b [Z from dedx in I Calo plane vs beta].. ............ %f \n",ZCalo_dedx_b);
  printf(" ZCalo_dedx_defl [Z from dedx in I Calo plane vs deflection....... %f \n",ZCalo_dedx_defl);
  printf(" ZCalo_Nmin1_defl [Z from truncated mean (N-1 pl) vs deflection].. %f \n",ZCalo_Nmin1_defl);
  printf("========================================================================\n");
  //
};

void CaloNuclei::Delete(){
  Clear();
  //delete this;
};


void CaloNuclei::Process(){
  Process(0);
};

void CaloNuclei::Process(Int_t ntr){
  //  
  if ( !L2 ){
    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
    printf(" ERROR: CaloNuclei variables not filled \n");
    return;
  };
  //
  Bool_t newentry = false;
  //
  if ( L2->IsORB() ){
    if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || ntr != sntr ){
      newentry = true;
      OBT = L2->GetOrbitalInfo()->OBT;
      PKT = L2->GetOrbitalInfo()->pkt_num;
      atime = L2->GetOrbitalInfo()->absTime;
      sntr = ntr;
    };
  } else {
    newentry = true;
  };
  //
  if ( !newentry ) return;
  //
  tr = ntr;
  //
  if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
  //
  Clear();
  //
  if ( debug ) printf(" Always calculate stdedx1 \n");
  //
  // Always calculate stdedx1 and maxrel
  //
  Int_t cont=0;
  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 mip = 0.;
  for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
    //
    mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
    //
    // put in vfpl vector the energy release on the first plane
    //
    if ( strip != -1 && view == 1 && plane == 0 ) {
      stfpl[indx] = strip;
      vfpl[indx] = mip;
      indx++;
    };
    //
  };
  //
  if ( debug ) printf(" find energy released along the strip of maximum on the first plane and on the two neighbour strips \n");
  //
  // 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<<stdedx1<<"      "<<maxrel<<"\n";
  //
  if ( debug ) printf(" if ( !usetrack ) return: usetrack %i ntr %i \n",usetrack,ntr);
  // 
  //
  //  if ( !usetrack ) return;
  //
  PamTrack *ptrack = 0;
  CaloTrkVar *track = 0;
  //
  if ( usetrack ){
    if ( ntr >= 0 ){
      ptrack = L2->GetTrack(ntr);
      if ( ptrack ) track = ptrack->GetCaloTrack();
    } else {
      track = L2->GetCaloStoredTrack(ntr);
    };
    //
    if ( !track && ntr >= 0 ){
      printf(" ERROR: cannot find any track!\n");
      printf(" ERROR: CaloNuclei variables not completely filled \n");
      return;   
    };
  } else {
    if ( ntr >= 0 ){
      if ( debug ) printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr);
      if ( debug ) printf(" ERROR: CaloNuclei variables not completely filled \n");
      return;       
    };
  };
  //
  //  Float_t defethr = 6. * 0.90;
  Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.;
  //
  // Calculate dedx1 and dedx3
  //
  for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
    //
    mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
    //
    if ( ntr >= 0 ){
      //
      if ( strip != -1 && 
	   view == 1   && 
	   plane == 0  &&
	   ( strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) || strip == (track->tibar[0][1]) ) 
	   && true ){      
	dedx1 += mip;
      };
      if ( strip != -1 && 
	   (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
	    ( view == 0 && plane == 0 )) &&
	   (( view == 0 && ( strip == track->tibar[0][0] || strip == (track->tibar[0][0]-1) || strip == (track->tibar[0][0]-2) )) ||
	    ( view == 1 && ( strip == track->tibar[0][1] || strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) )) ||
	    ( view == 1 && ( strip == track->tibar[1][1] || strip == (track->tibar[1][1]-1) || strip == (track->tibar[1][1]-2) ))) &&
	   true ){
	dedx3 += mip;
      };
    } else {
      //
      if ( strip != -1 && 
	   view == 1   && 
	   plane == 0  &&
	   ( strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) || strip == (L2->GetCaloLevel2()->cibar[0][1]) ) 
	   && true ){      
	dedx1 += mip;
      };
      if ( strip != -1 && 
	   (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
	    ( view == 0 && plane == 0 )) &&
	   (( view == 0 && ( strip == L2->GetCaloLevel2()->cibar[0][0] || strip == (L2->GetCaloLevel2()->cibar[0][0]-1) || strip == (L2->GetCaloLevel2()->cibar[0][0]-2) )) ||
	    ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[0][1] || strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) )) ||
	    ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[1][1] || strip == (L2->GetCaloLevel2()->cibar[1][1]-1) || strip == (L2->GetCaloLevel2()->cibar[1][1]-2) ))) &&
	   true ){
	dedx3 += mip;
      };
    };
    //    
  };
  //
  dedx3 /= 3.;
  //  Float_t mesethr = dedx1 * 0.90;
  Float_t mesethr = 0.;
  if ( dedx1 > 0. ) mesethr = (sqrt(dedx1) - 0.50)*(sqrt(dedx1) - 0.50);
  Bool_t aldone = false;
  //
 retry:
  // 
  if ( debug ) printf("retry\n");
  //
  interplane = 0;
  //
  ethr = TMath::Max(defethr,mesethr);
  //
  // Find the interaction plane "interplane"
  //
  Int_t gapth = 3;
  Int_t nhit[2] = {0,0};
  Int_t splane[2] = {-1,-1};
  Int_t sview[2] = {-1,-1};
  Int_t interpl[2] = {-1,-1};
  Int_t interv[2] = {-1,-1};
  Bool_t wmulthit[2] = {false,false};
  Bool_t wgap[2] = {false,false};
  Int_t ii = 0;
  while ( ii<L2->GetCaloLevel1()->istrip ){
    //
    mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
    //
    if ( ntr >= 0 ){
      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
	   ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ) 
	   && true ){      
	if ( debug ) printf(" inside loop: ii %i mip %f view %i plane %i strip %i tibar %i nhit %i splane %i sview %i \n",ii,mip,view,plane,strip,track->tibar[plane][view]-1,nhit[view],splane[view],sview[view]);
	interpl[view] = plane;
	interv[view] = view;
	if ( splane[view] != plane || sview[view] != view ){
	  if ( nhit[view] > 1 ){
	    wmulthit[view] = true;
	    //	  if ( splane[view] == -1 ) splane[view] = 0; //
	    //	  if ( sview[view] == -1 ) sview[view] = view; //
	    interpl[view] = splane[view];
	    interv[view] = sview[view];	  
	};
	  if ( plane > splane[view]+gapth ){
	    wgap[view] = true;
	    //	  if ( splane[view] == -1 ) splane[view] = 0;//
	    //	  if ( sview[view] == -1 ) sview[view] = view; //
	    interpl[view] = splane[view];
	    interv[view] = sview[view];
	  };
	  splane[view] = plane;
	  sview[view] = view;
	  nhit[view] = 1;
	} else {
	  nhit[view]++;
	};
      };
    } else {
      if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
	   ( strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ) 
	   && true ){      
	if ( debug ) printf(" inside loop: ii %i mip %f view %i plane %i strip %i cibar %i nhit %i splane %i sview %i \n",ii,mip,view,plane,strip,L2->GetCaloLevel2()->cibar[plane][view]-1,nhit[view],splane[view],sview[view]);
	interpl[view] = plane;
	interv[view] = view;
	if ( splane[view] != plane || sview[view] != view ){
	  if ( nhit[view] > 1 ){
	    wmulthit[view] = true;
	    //	  if ( splane[view] == -1 ) splane[view] = 0; //
	    //	  if ( sview[view] == -1 ) sview[view] = view; //
	    interpl[view] = splane[view];
	    interv[view] = sview[view];	  
	};
	  if ( plane > splane[view]+gapth ){
	    wgap[view] = true;
	    //	  if ( splane[view] == -1 ) splane[view] = 0;//
	    //	  if ( sview[view] == -1 ) sview[view] = view; //
	    interpl[view] = splane[view];
	    interv[view] = sview[view];
	  };
	  splane[view] = plane;
	  sview[view] = view;
	  nhit[view] = 1;
	} else {
	  nhit[view]++;
	};
      };
    };
    //
    ii++;
    //    
  };
  //
  if (debug ) printf("conversion interpl %i interv %i multhit %i interplane %i \n",interpl[0],interv[0],multhit,interplane);
  Int_t winterplane[2] = {-1,-1};
  //
  for ( Int_t view = 0; view < 2; view++){
    //
    if ( nhit[view] > 1 && !wmulthit[view] && !wgap[view] ){
      wmulthit[view] = true;
      interpl[view] = splane[view];
      interv[view] = sview[view];
    };
    //
    if ( wmulthit[view] ) multhit = true;
    if ( wgap[view] ) gap = true;
    //
    // convert view and plane number of interaction plane into number of available dE/dx measurements before the interaction plane
    //    
    if ( interpl[view] >= 0 ) {
      if ( interv[view] == 0 ){
	winterplane[view] = (1 + interpl[view]) * 2;
      } else {
	winterplane[view] = (1 + interpl[view]) + (1 + interpl[view] - 1);
      };
      if ( wmulthit[view] ) winterplane[view]--;
    };
  };
  if ( winterplane[0] > 0 && winterplane[1] > 0 ){
    if ( multhit ){
      interplane = TMath::Min(winterplane[0],winterplane[1]);
    } else {
      interplane = TMath::Max(winterplane[0],winterplane[1]);
    };
  } else {
    if ( !winterplane[0] || !winterplane[1] ){
      interplane = 0;
    } else {
      interplane = TMath::Max(winterplane[0],winterplane[1]);
    };
  };
  //  
  if ( debug ) printf("2conversion interpl %i interv %i multhit %i interplane %i \n",interpl[1],interv[1],multhit,interplane);
  if ( debug ) printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);
  //
  Int_t ipl = 0;
  if ( interplane > 0 ){
    //
    // Calculate preq, postq, qpremean
    //
    cont++;
    ii = 0;
    Int_t ind = -1;
    Int_t qsplane = -1;
    Int_t qsview = -1;
    Int_t ind2 = -1;
    Int_t qsplane2 = -1;
    Int_t qsview2 = -1;
    Float_t qme[200];
    memset(qme,0,200*sizeof(Float_t));
    Float_t qme2[2112];
    memset(qme2,0,2112*sizeof(Float_t));
    //
    while ( ii<L2->GetCaloLevel1()->istrip ){
      //
      mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip);    
      //
      if ( strip != -1 ){
	if ( view == 0 ){
	  ipl = (1 + plane) * 2;
	} else {
	  ipl = (1 + plane) + (1 + plane - 1 );
	};
	if ( ipl > interplane ){
	  postq += mip;
	} else {
	  preq += mip;
	  if ( ntr >= 0 ){
	    if (  strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ){
	      if ( qsplane != plane || qsview != view ){
		qsplane = plane;
		qsview = view;
		ind++;
		if ( debug && ind > 199 ) printf(" AAAGH!! \n");
		qme[ind] = 0.;
	      };
	      qme[ind] += mip;
	    };
	    for ( Int_t ns = 0; ns < R ; ns++){
	      Int_t ms =  track->tibar[plane][view] - 1 - ns + (R - 1)/2;
	      if ( strip == ms ){
		if ( qsplane2 != plane || qsview2 != view ){
		  qsplane2 = plane;
		  qsview2 = view;
		  ind2++;
		  if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
		  qme2[ind2] = 0.;
		};
		qme2[ind2] += mip;
	      };
	    };	
	  } else {
	    if (  strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ){
	      if ( qsplane != plane || qsview != view ){
		qsplane = plane;
		qsview = view;
		ind++;
		if ( debug && ind > 199 ) printf(" AAAGH!! \n");
		qme[ind] = 0.;
	      };
	      qme[ind] += mip;
	    };
	    for ( Int_t ns = 0; ns < R ; ns++){
	      Int_t ms =  L2->GetCaloLevel2()->cibar[plane][view] - 1 - ns + (R - 1)/2;
	      if ( strip == ms ){
		if ( qsplane2 != plane || qsview2 != view ){
		  qsplane2 = plane;
		  qsview2 = view;
		  ind2++;
		  if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n");
		  qme2[ind2] = 0.;
		};
		qme2[ind2] += mip;
	      };
	    };	
	  };
	};      
	//
      };
      //
      ii++;
      //    
    };


 //
    // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean...
    //    
    if ( debug ){
      for (Int_t l=0; l < interplane; l++){
	printf(" qme[%i] = %f \n",l,qme[l]);  
      };
    };
    //
    Long64_t work[200];
    ind = 0;
    Int_t l = 0;
    Int_t RN = 0;
    Float_t qm = 0.;
    Float_t qm2 = 0.;
    //
    Float_t qmt = ethr*0.8; // *0.9
    //
    Int_t uplim = TMath::Max(3,N);
    Int_t uplim2 = interplane-1;
    //
    while ( l < uplim && ind < interplane ){
      qm = TMath::KOrdStat((Long64_t)interplane,qme,(Long64_t)ind,work);
      if ( qm >= qmt ){	
	if ( l < 3 ){
	  qpremean += qm;
	  RN++;
	};
	l++;
	if ( debug ) printf(" value no %i qm %f qmt %f \n",l,qm,qmt); 
      };
      ind++;
    };
    //
    qpremean /= (Float_t)RN;
    ind = 0;
    l = 0;
    RN = 0;
    while ( l < uplim && ind < interplane ){
      qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
      if ( qm2 >= qmt ){	
	if ( l < N ){
	  qpremeanN += qm2;
	  RN++;
	};
	l++;
	if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN); 
      };
      ind++;
    };
    ////////////////////////////////////
    //to calculate qNmin1///////////////
    ///////////////////////////////////
    //values under threshold
    qm2=0;
    ind = 0;
    l = 0;
    RN = 0;
    S2=0;
    while ( l < uplim2 && ind<interplane){
      qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
      if ( qm2 < qmt ) S2++;
      ind++;
    }
    qm2=0;
    ind = 0;
    l = 0;
    RN = 0;
    while ( l < uplim2 && ind < interplane ){
      qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work);
      if ( qm2 >= qmt ){	
	if ( l < (interplane - 1 - S2)){  
	  qNmin1_w += qm2;
	  RN++;
	};
	l++;
	if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN); 
      };
      ind++;
    };
    qpremeanN /= (Float_t)RN;
    qNmin1_w /= (Float_t)RN;
    UN = RN;
    ///////set qNmin1 definition///////////
    if (interplane==1 || interplane==2){  
      if (dedx1>0) qNmin1=dedx1;
      else if (stdedx1>0) qNmin1=stdedx1;	       
    }
    else if (interplane > 2){
      qNmin1 =  qNmin1_w;
    }
    ////////////////////////////////////
    //////////////////////////////////
    //
    if ( debug ) printf(" charge is %f \n",sqrt(qpremean));
    //
    if ( mesethr != ethr && interplane >= 3 && !aldone ){
      Float_t mesethr2 = (sqrt(qpremean) - 0.50)*(sqrt(qpremean) - 0.50);
      if ( mesethr2 < mesethr*0.90 ){
	mesethr = (sqrt(dedx1) - 0.25)*(sqrt(dedx1) - 0.25);
      } else {
	mesethr = mesethr2;
      };
      aldone = true;
      if ( mesethr > defethr ){
	interplane = 0;
	preq = 0.;
	postq = 0.;
	qpremean = 0.; 
	qpremeanN = 0.; 
	qNmin1 = 0;
	multhit = false;
	gap = false;
	goto retry;
      };
    };
  };



//=======================================================================
//===========    charge determination stdedx1 vs. beta    ===============
//======================     Siegen method    ===========================
//=======================================================================

// 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;

Float_t charge = 1000.;
Float_t beta = 100.;

//------- First try track dependent beta
    if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
    PamTrack *TRKtrack = L2->GetTrack(0);
    if (fabs(TRKtrack->GetToFTrack()->beta[12]) < 100.) beta = fabs(TRKtrack->GetToFTrack()->beta[12]);
                                                }
//------- If no beta found, try standalone beta
    if (beta == 100.) {
    ToFTrkVar *ttrack = L2->GetToFStoredTrack(-1);
    beta = fabs(ttrack->beta[12]);
                      }

    if (beta<2.) {  // it makes no sense to allow beta=5 or so...

    Float_t  mip=0;
    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 < 100


    charge_siegen1 = charge;

//=======================  end charge Siegen  ===========================


// //=======================================================================
// //===========    charge determination Maximum release  vs. beta    ===============
// //======================     Rome method    ===========================
// //=======================================================================
    
    Float_t D0[7] = {0, 3 , 4 , 5 , 6, 8, 90};
    Float_t E1[7] = {0 ,923.553 , 659.842, 1113.97, 3037.25, 3034.84, 0};
    Float_t E2[7] = {0 ,6.92574 , 5.08865, 5.29349, 6.41442, 5.52969, 0};
    Float_t E3[7] = {0 ,9.7227  , 13.18, 23.5444, 38.2057, 63.6784, 80000};

    Float_t xx1[7],yy1[7];
    n1 = 7;
    
    charge = 1000.;
    mip=0;
    

    if (beta<2.) {  // it makes no sense to allow beta=5 or so...
      
      
      mip=maxrel;

      if (mip>0)   {
	Float_t ym = mip;
	Float_t xb = beta;
	
	for ( Int_t jj=0; jj<n1; jj++ ){
	  xx1[jj] = E1[jj]*exp(-E2[jj]*xb)+E3[jj];
	  yy1[jj] = D0[jj]*D0[jj] ;
	}
	
	TGraph *gr1 = new TGraph(n1,xx1,yy1);
	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 < 100


    ZCalo_maxrel_b = charge;

    //=======================  end charge Rome: maxril vs beta  ===========================



// =======================================================================
// ===========    charge determination dedx  vs. beta    ===============
// ======================     Rome method    ===========================
// =======================================================================
    
    Float_t F0[7] = {0.,3. ,4., 5. , 6., 8, 90};
    Float_t G1[7] = {0 ,642.935 , 848.684, 1346.05, 3238.82, 3468.6, 0};
    Float_t G2[7] = {0 ,6.2038 , 5.51723, 5.65265, 6.54089, 5.72723, 0};
    Float_t G3[7] = {0 ,9.2421 , 13.9858, 25.3912, 39.6332, 64.5674, 80000};

  
    charge = 1000.;
    mip=0;


    if (beta<2.) { // it makes no sense to allow beta=5 or so...
      
    
      if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
	mip=dedx1;
      } 
      if (mip==0) mip=stdedx1;
      

      if (mip>0)   {
	
	Float_t ym = mip;
	Float_t xb = beta;
	
	for ( Int_t jj=0; jj<n1; jj++ ){
 	 xx1[jj] = G1[jj]*exp(-G2[jj]*xb)+G3[jj];
	 yy1[jj] = F0[jj]*F0[jj] ;
	}
	
	TGraph *gr1 = new TGraph(n1,xx1,yy1);
	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 < 100
    
    ZCalo_dedx_b = charge;

    //=======================  end charge Rome: dedx vs beta  ===========================
    
    
//=======================================================================
//===========    charge determination dedx  vs. defl    ===============
//======================     Rome method    ===========================
//=======================================================================
    
    //new
     Float_t H0[7] = {0, 3 , 4 , 5 , 6, 8, 90 };
     Float_t I1[7] = {0 , 56.1019, 101.673, 109.225, 150.599, 388.531, 0};
     Float_t I2[7] = {0 , -12.5581, -22.5543, -15.9823, -28.2207, -93.6871, 0};
     Float_t I3[7] = {0 , 11.6218, 19.664, 32.1817, 45.7527, 84.5992, 80000};
    

//     Float_t H0[7] = {0, 3 , 4 , 5 , 6, 8, 90 };
//     Float_t I1[7] = {0 , 56.1019, 101.673, 155, 150.599, 388.531, 0};
//     Float_t I2[7] = {0 , -12.5581, -22.5543, -35.6217, -28.2207, -93.6871, 0};
//     Float_t I3[7] = {0 , 11.6218, 19.664, 34.3311, 45.7527, 84.5992, 8000};


    
    charge = 1000.;
    mip=0;
    Float_t defl=0;
    

    if (beta<2.) { // it makes no sense to allow beta=5 or so...

      if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
	PamTrack *TRKtrack = L2->GetTrack(0);
	mip=dedx1;
	if (mip==0) mip=stdedx1;
	defl=TRKtrack->GetTrkTrack()->al[4];
	
	
	if (mip>0 && defl<0.7 && defl>0)   {
	  
	  Float_t ym = mip;
	  Float_t xb = defl;
	  
	  for ( Int_t jj=0; jj<n1; jj++ ){
	    xx1[jj] = I1[jj]*xb*xb+I2[jj]*xb+I3[jj];
	    yy1[jj] = H0[jj]*H0[jj] ;
	  }
	  
	  TGraph *gr1 = new TGraph(n1,xx1,yy1);
	  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 && defl<0.5 && defl>0)
      }//Ntrack>=1
    } //beta < 100

    ZCalo_dedx_defl = charge;
 
//=======================  end charge Rome: dedx vs defl  ===========================


//============================================================================================
//===========    charge determination Truncated mean (N-1 planes) vs. defl ===================
//================================     Rome method    ========================================
//============================================================================================
    
    Float_t L0[7] = {0, 3 , 4 , 5 , 6, 8, 90};
    Float_t M1[7] = {0 , 63.0145, 120.504, 173.663, 245.33, 236.517, 0};
    Float_t M2[7] = {0 , -15.005, -31.0635, -39.4988, -60.5011, -46.3992, 0};
    Float_t M3[7] = {0 , 12.5037, 22.8652, 35.2907, 51.4678, 86.4155, 8000}; 

    charge = 1000.;
    mip=0;
   

    if (beta<2.) { // it makes no sense to allow beta=5 or so...

      if( L2->GetTrkLevel2()->GetNTracks()>=1 ){
	mip=qNmin1;
	
	if (mip>0 && defl<0.7 && defl>0)   {
	  
	  Float_t ym = mip;
	  Float_t xb = defl;
	  
	  for ( Int_t jj=0; jj<n1; jj++ ){
	    xx1[jj] = M1[jj]*xb*xb+M2[jj]*xb+M3[jj];
	    yy1[jj] = L0[jj]*L0[jj] ;
	  }
	  
	  TGraph *gr1 = new TGraph(n1,xx1,yy1);
	  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 && defl<0.5 && defl>0)
      }//Ntrack>=1
    } //beta < 100
    
    ZCalo_Nmin1_defl = charge;
 
//=======================  end charge Rome: Nmin1 vs defl  ===========================





  //
  if ( debug ) this->Print();
  if ( debug ) printf(" esci \n");
  //
};