#include <CaloNuclei.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;
  //
};

void CaloNuclei::Clear(){
  //
  interplane = 0;
  preq = 0.;
  postq = 0.;
  dedx1 = 0.;
  dedx3 = 0.;
  qpremean = 0.; 
  qpremeanN = 0.; 
  //
  multhit = false;
  gap = false;
  //
};

void CaloNuclei::Print(){
  //
  printf("==============================\n");
  printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
  printf(" interplane: %i\n",interplane);
  printf(" ethr: %f \n",ethr); 
  printf(" dedx1: %f \n",dedx1);
  printf(" dedx3: %f \n",dedx3);
  printf(" multhit: %i \n",multhit);
  printf(" gap: %i \n",gap);
  printf(" preq: %f \n",preq);
  printf(" postq: %f \n",postq);
  printf(" qpremean: %f \n",qpremean); 
  printf(" N: %i \n",N); 
  printf(" qpremeanN: %f \n",qpremeanN); 
  printf("==============================\n");
  //
};

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


void CaloNuclei::Process(){
  //
  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 ){
      newentry = true;
      OBT = L2->GetOrbitalInfo()->OBT;
      PKT = L2->GetOrbitalInfo()->pkt_num;
      atime = L2->GetOrbitalInfo()->absTime;
    };
  } else {
    newentry = true;
  };
  //
  if ( !newentry ) return;
  //
  //  printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
  //
  Clear();
  //
  PamTrack *track = 0;
  track = L2->GetTrack(0);
  //
  Int_t view = 0;
  Int_t plane = 0;
  Int_t strip = 0;
  Float_t mip = 0.;
  //  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 ( strip != -1 && 
	 view == 1   && 
	 plane == 0  &&
	 ( strip == (track->GetCaloTrack()->tibar[0][1]-1) || strip == (track->GetCaloTrack()->tibar[0][1]-2) || strip == (track->GetCaloTrack()->tibar[0][1]) ) 
	 && true ){      
      dedx1 += mip;
    };
    if ( strip != -1 && 
	 (( view == 1 && ( plane == 0 || plane == 1 ) ) ||
	 ( view == 0 && plane == 0 )) &&
	 (( view == 0 && ( strip == track->GetCaloTrack()->tibar[0][0] || strip == (track->GetCaloTrack()->tibar[0][0]-1) || strip == (track->GetCaloTrack()->tibar[0][0]-2) )) ||
	  ( view == 1 && ( strip == track->GetCaloTrack()->tibar[0][1] || strip == (track->GetCaloTrack()->tibar[0][1]-1) || strip == (track->GetCaloTrack()->tibar[0][1]-2) )) ||
	  ( view == 1 && ( strip == track->GetCaloTrack()->tibar[1][1] || strip == (track->GetCaloTrack()->tibar[1][1]-1) || strip == (track->GetCaloTrack()->tibar[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:
  // 
  //  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 ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] &&
	 ( strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) ) 
	 && true ){      
      //      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->GetCaloTrack()->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;
	  interpl[view] = splane[view];
	  interv[view] = sview[view];
	  // 	  ii = L2->GetCaloLevel1()->istrip;
	};
	if ( plane > splane[view]+gapth ){
	  wgap[view] = true;
	  interpl[view] = splane[view];
	  interv[view] = sview[view];
	  //	  ii = L2->GetCaloLevel1()->istrip;
	};
	splane[view] = plane;
	sview[view] = view;
	nhit[view] = 1;
      } else {
	nhit[view]++;
      };
    };
    //
    ii++;
    //    
  };
  //
  //  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]);
    };
  };
  //  
  //  printf("2conversion interpl %i interv %i multhit %i interplane %i \n",interpl[1],interv[1],multhit,interplane);
  //  printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]);
  //
  Int_t ipl = 0;
  if ( interplane > 0 ){
    //
    // Calculate preq, postq, qpremean
    //
    ii = 0;
    Int_t ind = -1;
    Int_t qsplane = -1;
    Int_t qsview = -1;
    Float_t qme[200];
    memset(qme,0,200*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 (  strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) ){
	    if ( qsplane != plane || qsview != view ){
	      qsplane = plane;
	      qsview = view;
	      ind++;
	      if ( ind > 199 ) printf(" AAAGH!! \n");
	      qme[ind] = 0.;
	    };
	    qme[ind] += mip;
	  };
	};      
	//
      };
      //
      ii++;
      //    
    };
    //
    // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean...
    //    
    //    for (Int_t l=0; l < interplane; l++){
    //       printf(" qme[%i] = %f \n",l,qme[l]);  
    //       //      if ( qme[ind] < 6.25 ) qme[ind] = 10000.;      
    //     };
    Long64_t work[200];
    ind = 0;
    Int_t l = 0;
    Float_t qm = 0.;
    Int_t uplim = TMath::Max(3,N);
    while ( l < uplim && ind < interplane ){
      qm = TMath::KOrdStat(interplane,qme,ind,work);
      if ( qm >= mesethr*0.9 ){	
	if ( l < 3 ) qpremean += qm;
	if ( l < N ) qpremeanN += qm;
	l++;
	//	printf(" value no %i qm %f \n",l,qm); 
      };
      ind++;
    };
    //
    qpremean /= l;
    qpremeanN /= N;
    //
    //  printf(" charge is %f \n",sqrt(qpremean));
    //
    //    printf("preq postq\n");    
    //
    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 ) goto retry;
    };
  };
  //
  //  printf(" esci \n");
  //
};