#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;
  R = 3;
  //
  debug = false;
  usetrack = true;
  //
};

void CaloNuclei::Clear(){
  //
  tr = 0;
  sntr = 0;
  interplane = 0;
  preq = 0.;
  postq = 0.;
  dedx1 = 0.;
  dedx3 = 0.;
  qpremean = 0.; 
  qpremeanN = 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("========================================================================\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
  //
  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,stfpl);
    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];
    };
  } else {
    stdedx1 = 0.;
  };
  //
  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 ){
      printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr);
      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
    //
    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);
    //
    while ( l < uplim && ind < interplane ){
      qm = TMath::KOrdStat(interplane,qme,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(interplane,qme2,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++;
    };
    //
    qpremeanN /= (Float_t)RN;
    //
    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.; 
	multhit = false;
	gap = false;
	goto retry;
      };
    };
  };
  //
  if ( debug ) this->Print();
  if ( debug ) printf(" esci \n");
  //
};