//*-- AUTHOR : Y.C.Pachmayer
//*-- Modified :
///////////////////////////////////////////////////////////////////////////////
//
// HMdcCosmicsCandidate
//
//
///////////////////////////////////////////////////////////////////////////////


using namespace std;
#include "hmdccosmicscandidate.h"
#include "TH1.h"
#include "TH2.h"
#include "TROOT.h"
#include "TMath.h"
#include "hades.h"
#include "hevent.h"
#include "hruntimedb.h"
#include "haddef.h"
#include "hiterator.h"
#include "hmdcdef.h"
#include "hmdctrackddef.h"
#include "hmdcdetector.h"
#include "tofdef.h"
#include "hspectrometer.h"
#include "htofgeompar.h"
#include "hmdccal1sim.h"
#include "htofhit.h"
#include "htofcluster.h"
#include "hmdcclussim.h"
#include "hmdcsizescells.h"
#include "hcategory.h"
#include "hgeomtransform.h"
#include "hmdcgetcontainers.h"
#include <iostream>
#include <fstream>

ClassImp(HMdcCosmicsCandidate)
ClassImp(HCosmicCalibEvSkip)


HMdcCosmicsCandidate::HMdcCosmicsCandidate(const Text_t *name,const Text_t *title) 
    : HReconstructor(name,title) {
  initVariables();
}

HMdcCosmicsCandidate::HMdcCosmicsCandidate(void)
    : HReconstructor("HMdcCosmicsCandidate","HMdcCosmicsCandidate") {
  initVariables();
}

HMdcCosmicsCandidate::~HMdcCosmicsCandidate(void) {
  if (tofIter) delete tofIter;
  tofIter  = NULL;
  HMdcSizesCells::deleteCont();
  HMdcGetContainers::deleteCont();
  if(histFile != NULL) delete histFile;
  histFile = NULL;
}

void HMdcCosmicsCandidate::initVariables(void) {
  returnFlag     = kSkipEvent;
  pCatTofHit     = NULL;
  pCatTofCluster = NULL;
  tofIter        = NULL;
  pCatMdcCal1    = NULL;
  hsBeta         = NULL;
  plStat         = NULL;
  for(Int_t s=0;s<6;s++) for(Int_t m=0;m<4;m++) {
    hsDCell[s][m] = NULL;
    for(Int_t l=0;l<6;l++) hsDCellL[s][m][l] = NULL;
  }

  makeHists      = kFALSE;
  histFile       = NULL;
  histFileName   = "./cosmCandHists.root";
  histFileOption = "NEW";
  setParContStat = kFALSE;
  zTofShift      = 0.;
  
  // Cuts:
  betaMin        = 0.5;
  betaMax        = 1.5;
  for(Int_t m=0; m<4; m++) dCellWind[m]    = 15.;
  for(Int_t s=0;s<6;s++)   nModInSecCut[s] = 2; // Minimum 2 modules
  for(Int_t s=0;s<6;s++) for(Int_t m=0; m<4; m++) nLayInModCut[s][m] = 3;
  nEventsTot   = 0;
  nEventsSel   = 0;
  nEventsSelSc = 0;
  for(Int_t i=0;i<9;i++) {
    nSelTr[i]       = 0;
    trackScaling[i] = i<4 ? 0. : 1.1;
    nSelTrSc[i]     = 0;
  }
}

void HMdcCosmicsCandidate::setBetaCut(Float_t min,Float_t max) {
  if(max>0. && min>=0. && max>min) {
    betaMin = min;
    betaMax = max;
  } else {
    Error("setBetaCut","Error in beta cut setting: min=%f max=%f",min,max);
    exit(1);
  }
}

void HMdcCosmicsCandidate::setDCellWindow(Float_t wM1,Float_t wM2,Float_t wM3,Float_t wM4) {
  if(wM1>=0.) dCellWind[0] = wM1;
  if(wM2>=0.) dCellWind[1] = wM2;
  if(wM3>=0.) dCellWind[2] = wM3;
  if(wM4>=0.) dCellWind[3] = wM4;
}

void HMdcCosmicsCandidate::setNLayInModCut(Int_t sec,Int_t lm1,Int_t lm2,Int_t lm3,Int_t lm4) {
  if(sec<0 || sec>5) return;
  nLayInModCut[sec][0] = lm1;
  nLayInModCut[sec][1] = lm2;
  nLayInModCut[sec][2] = lm3;
  nLayInModCut[sec][3] = lm4;
}
  
Bool_t HMdcCosmicsCandidate::init(void) {
  cout << " HMdcCosmicsCandidate " << endl;
  if(!gHades) return kFALSE;

  sizes = HMdcSizesCells::getObject();
  pGetCont = HMdcGetContainers::getObject();

  pCatTofCluster = gHades->getCurrentEvent()->getCategory(catTofCluster);
  if(pCatTofCluster != NULL) {
    tofIter = (HIterator *)pCatTofCluster->MakeIterator("native");
    if(!tofIter) return kFALSE;
  } else {
    pCatTofHit = gHades->getCurrentEvent()->getCategory(catTofHit);
    if (!pCatTofHit) {
      Error("HMdcCosmicsCandidate::init()","NO HTofHit found!");
      return kFALSE;
    }
    tofIter = (HIterator *)pCatTofHit->MakeIterator("native");
    if(!tofIter) return kFALSE;
  }
    
  pCatMdcCal1 = gHades->getCurrentEvent()->getCategory(catMdcCal1);
  if (pCatMdcCal1 == NULL) {
    Error("HMdcCosmics::init()","NO HMdcCal1 found!");
    return kFALSE;
  }
  pCatMdcClus = pGetCont->getCatMdcClus(kTRUE);
  if(pCatMdcClus == NULL) {
    Error("HMdcCosmics::init()","NO HMdcClus found!");
    return kFALSE;
  }

  //  MDC setup:
  HMdcDetector*  fMdcDet = pGetCont->getMdcDetector();
  for(Int_t s=0; s<6; s++) for(Int_t m=0; m<4; m++) {
    if(fMdcDet->getModule(s,m) == 0) mdcsetup[s][m] = 0;
    else                             mdcsetup[s][m] = 1;
  }

  printParam();
  if(makeHists && !createHists()) return kFALSE;
  return kTRUE;

}

Bool_t HMdcCosmicsCandidate::reinit(void) { 
  if(!sizes->initContainer()) return kFALSE;
  return kTRUE;
}

void HMdcCosmicsCandidate::printParam(void) {
  printf("------------------------------------------------------------\n");
  printf("  HMdcCosmicsCandidate:\n  Tof");
  if(pCatTofCluster != NULL)  printf("Clusters");
  else if(pCatTofHit != NULL) printf("Hits");
  printf(" are used for cosmic track finder.\n");
  printf("  Beta cut                   : betaMin=%g, betaMax=%g\n",betaMin,betaMax);
  if(zTofShift != 0.f) printf("  Shift of TofHits Z position: %.1f mm\n",zTofShift);
  printf("  Cells cut (MDC I,II,III,IV): %.1f, %.1f, %.1f, %.1f\n",
         dCellWind[0],dCellWind[1],dCellWind[2],dCellWind[3]);
  if(makeHists) printf("  (See histograms in the file)\n");
  printf("  Filtering:\n");
  printf("  Minimal number of hitted modules in sector[1-6]:");
  for(Int_t s=0;s<6;s++) printf(" %2i",nModInSecCut[s]);
  printf("\n");
  for(Int_t s=0;s<6;s++) {
    printf("  Sector %i. Min.number of wires in module[I-IV]:",s+1);
    for(Int_t m=0; m<4; m++) printf(" %2i",nLayInModCut[s][m]);
    printf("\n");
  }
  printf("------------------------------------------------------------\n");
}

void HMdcCosmicsCandidate::setNModInSecCut(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5,Int_t s6) {
  nModInSecCut[0] = s1;
  nModInSecCut[1] = s2;
  nModInSecCut[2] = s3;
  nModInSecCut[3] = s4;
  nModInSecCut[4] = s5;
  nModInSecCut[5] = s6;
}

Int_t HMdcCosmicsCandidate::execute(void) {
  
  if(nEventsTot==0 && setParContStat) {
    gHades->getRuntimeDb()->setContainersStatic();
    gHades->setEnableCloseInput(kFALSE);
  }
  
  nEventsTot++;

  Bool_t isGeant = HMdcGetContainers::isGeant();
  
  tofMult = 0;
  if(pCatTofCluster != NULL) {
    HTofCluster* tofCluster;
    tofIter->Reset();
    while ( (tofCluster=(HTofCluster *)tofIter->Next()) != 0 ) {
      Int_t tofSector = tofCluster->getSector();
      if(&((*sizes)[tofSector]) == 0) return returnFlag; // No mdc in this sector
      nsec[tofMult] = tofSector;
      tof[tofMult]  = tofCluster->getTof();
      tofCluster->getXYZLab(xlab[tofMult],ylab[tofMult],zlab[tofMult]);
//if(zlab[tofMult]>490.) continue;
      zlab[tofMult] += zTofShift;

      tofMult++;
    }
  } else {
    HTofHit* tofHit;
    tofIter->Reset();
    while ( (tofHit=(HTofHit *)tofIter->Next()) != 0 ) {
      Int_t tofSector = tofHit->getSector();
      if(&((*sizes)[tofSector]) == 0) return returnFlag; // No mdc in this sector
      nsec[tofMult] = tofSector;
      tof[tofMult]  = tofHit->getTof();
      tofHit->getXYZLab(xlab[tofMult],ylab[tofMult],zlab[tofMult]);
      zlab[tofMult] += zTofShift;

      tofMult++;
    }
  }

  if (tofMult!=2 || nsec[0]==nsec[1]) return returnFlag;
  // Cut for track beta:
  Float_t dx = xlab[0]-xlab[1];
  Float_t dy = ylab[0]-ylab[1];
  Float_t dz = zlab[0]-zlab[1];
  Float_t path = sqrt(dx*dx + dy*dy + dz*dz);
  Float_t beta = path/TMath::Abs(tof[1]-tof[0])*1.e+6/TMath::C();
  if(makeHists) hsBeta->Fill(beta);
  if(beta < betaMin || beta > betaMax) return returnFlag;
  if(makeHists) {
    if(tof[0]<=tof[1]) plStat->Fill(nsec[0]+0.5,nsec[1]+0.5);
    else               plStat->Fill(nsec[1]+0.5,nsec[2]+0.5);
  }

//if(xlab[0]<0. && xlab[1]<0.) return returnFlag;
  
  //!!!!!!!!!! +++
  HMdcClus* prevClus = 0;
  Int_t     prevInd  = -1;

  Int_t nLayers      = 0;
  Int_t nModules     = 0;
  Int_t nWires       = 0;
  Int_t nModInSec[2] = {0,0}; // Number of modules in sector
          
  for(Int_t indTof=0;indTof<2;indTof++) {
    Int_t sec = nsec[indTof];
    HMdcSizesCellsSec&  sizessec = (*sizes)[sec];
    Double_t tofy1,tofz1;
    Double_t tofx1 = transTofHitToSec(0,sec, tofy1,tofz1);
    Double_t tofy2,tofz2;
    Double_t tofx2 = transTofHitToSec(1,sec, tofy2,tofz2);
    for(Int_t seg=0; seg<2;seg++) {
      //	if(seg==0) continue; // only outer chambers
//!!!		if(seg==1) continue; // only inner chambers

      HMdcList12GroupCells listCells;

      for(Int_t mod=seg*2; mod<(seg+1)*2;mod++) if(mdcsetup[sec][mod]>0) {
	HMdcSizesCellsMod& sizesmod = sizessec[mod];
	for(Int_t lay=0;lay<6;lay++) {
	  HMdcSizesCellsLayer& sizesLayer = sizesmod[lay];
          Float_t cell1, cell2;
	  sizesLayer.calcCrossedCells(tofx1, tofy1, tofz1,tofx2, tofy2, tofz2,cell1, cell2);
          Float_t cell12   = (cell1+cell2)/2;
          Short_t nCells   = sizesLayer.getNCells();

	  Int_t nFirstCell = (Int_t)(cell1 - dCellWind[mod]);
 	  Int_t nLastCell  = (Int_t)(cell2 + dCellWind[mod]);
          if(nFirstCell < 0)      nFirstCell = 0;
          if(nLastCell >= nCells) nLastCell  = nCells-1;

          Int_t nFCell = (Int_t)(cell1 - 30.);
          Int_t nLCell = (Int_t)(cell2 + 30.);
          if(nFCell < 0)       nFCell = 0;
          if(nLCell >= nCells) nLCell = nCells-1;

          Int_t lCells[nLCell-nFCell+1];
          Int_t nCs = 0;
	  for(Int_t cell=nFCell; cell<=nLCell;cell++) {
	    locCal1.set(4,sec,mod,lay,cell);
	    HMdcCal1* cal1 = (HMdcCal1*)pCatMdcCal1->getObject(locCal1);
	    if(cal1 == NULL) continue;
            if(cell>=nFirstCell && cell<=nLastCell) {
              lCells[nCs] = cell;
              nCs++;
            }
            // Fill histograms:
            if(makeHists) {
              hsDCell[sec][mod]->Fill(cell-cell12);
              hsDCellL[sec][mod][lay]->Fill(cell-cell12);
            }
//                 if(makeHists) {
//                   if     (cell<=cell1) hsDCell[sec][mod]->Fill(cell-cell1);
//                   else if(cell>=cell2) hsDCell[sec][mod]->Fill(cell-cell2);
//                   else                 hsDCell[sec][mod]->Fill(0.);
//                 }
	  }
          
          if(nCs>0) {
            Int_t layInd = lay+(mod%2)*6;
            if(nCs == 1) listCells.setTime(layInd, lCells[0], 1);
            else if(nCs==2 && lCells[1]-lCells[0]==1) {
              listCells.setTime(layInd, lCells[0], 1);
              listCells.setTime(layInd, lCells[1], 1);
            } else {
              Int_t   cM    = 0;
              Float_t dCell = TMath::Abs(Float_t(lCells[0]) - cell12);
              for(Int_t nc=1;nc<nCs;nc++) {
                Float_t dC = TMath::Abs(Float_t(lCells[nc]) - cell12);
                if(dC<dCell) {
                  dCell = dC;
                  cM    = nc;
                }
              }
              listCells.setTime(layInd, lCells[cM], 1);
              if(cM>0     && lCells[cM]-lCells[cM-1]==1) listCells.setTime(layInd,lCells[cM-1],1);
              if(cM<nCs-1 && lCells[cM+1]-lCells[cM]==1) listCells.setTime(layInd,lCells[cM+1],1);
//               for(Int_t nc=0;nc<nCs;nc++) {
//                 listCells.setTime(layInd, lCells[nc], 1);
//               }
            }
          }
	}
        Int_t nLay = listCells.getNLayersMod(mod%2);
        if(nLay>=nLayInModCut[sec][mod]) {
          nLayers += nLay;
          nModules++;
          nModInSec[indTof]++;
        }
      }
      Int_t nwr = listCells.getNCells();
      if(nwr == 0) continue;
      nWires += nwr;

      locClus.set(2,sec,seg);

      HMdcClus* fClus=(HMdcClus*)pCatMdcClus->getNewSlot(locClus, &index);
      if(fClus==0) {
        Error("execute","S.%i No HMdcClus slot available!",sec+1);
        continue;
      }
      if(isGeant) fClus = (HMdcClus*)(new(fClus) HMdcClusSim(listCells));
      else        fClus = new(fClus) HMdcClus(listCells); 

      fClus->setOwnIndex(index);

      //!!!!!!!!!!!!!!!!!!!!+++
      if(prevClus)  {
        prevClus->setIndexChilds(index,index);
        fClus->setIndexParent(prevInd);
      }
      prevClus = fClus;
      prevInd  = index;

      fClus->setSec(sec);
      fClus->setIOSeg(seg);
      fClus->setMod(-1);

      fClus->setTarget(tofx1, 20, tofy1, 30, tofz1, 30); 
      fClus->setXY(tofx2, 20, tofy2, 30);
      fClus->setPrPlane(0., 0., tofz2); // If parA=parB=0 then parD=Z
      if(isGeant) ((HMdcClusSim*)fClus)->calcTrList();
      // Cleaning:
      fClus->setMinCl(0,0);
      fClus->setRealLevel(0);
      fClus->setOwnIndex(index);
      fClus->setNBins(0);
    }
  }
  // Event selection:
  if(nModInSec[0]<nModInSecCut[nsec[0]])  return returnFlag;  // Main cut
  if(nModInSec[1]<nModInSecCut[nsec[1]])  return returnFlag;  // Main cut
  if(nModules<2 || nLayers<5 || nWires<5) return returnFlag;  // If main cut to low
  nSelTr[nModules]++;
  nEventsSel++;

  if(trackScaling[nModules] <= 0.)                                       return returnFlag;
  if(trackScaling[nModules] < 1.&& rndm.Rndm() > trackScaling[nModules]) return returnFlag;
  nSelTrSc[nModules]++;
  nEventsSelSc++;
  
  return 0;
 }


Float_t HMdcCosmicsCandidate::transTofHitToSec(Int_t tofInd, Int_t slocal,
                                               Double_t& tofy,Double_t& tofz) {
  // transformation of the tof hit from lab coordinate to sector coordinate system
  HMdcSizesCellsSec& sizessec = (*sizes)[slocal];
  
  HGeomVector A;    
  A.setXYZ(xlab[tofInd],ylab[tofInd],zlab[tofInd]); // vector A is lab coordinate system
  const HGeomTransform* vectortrans = sizessec.getLabTrans();
  HGeomVector a;
  a = vectortrans->transTo(A);             // here the transformation from lab to sec coordinate sytem is done
  Double_t tofx = a.getX();
  tofy = a.getY();
  tofz = a.getZ();
  return tofx;
}
  
void HMdcCosmicsCandidate::setHistFile(const Char_t* dir,const Char_t* fileNm,const Char_t* opt) {
  // the hist output file is set
  if(dir != NULL)                       histFileName   = dir;
  else                                  histFileName   = "./";
  if( !histFileName.EndsWith("/") )     histFileName  += '/';
  if(fileNm != NULL)                    histFileName  += fileNm;
  else                                  histFileName  += "cosmCandHists";
  if(histFileName.EndsWith(".hld"))     histFileName.Remove(histFileName.Length()-4);
  if( !histFileName.EndsWith(".root") ) histFileName  += ".root";
  if(opt>0)                             histFileOption = opt;
}

Bool_t HMdcCosmicsCandidate::createHists(void) {
  histFile = new TFile(histFileName.Data(),histFileOption.Data());
  if(histFile==NULL || !histFile->IsOpen()) {
    Error("createHists","Can't open file %s (opt=%s)!",histFileName.Data(),histFileOption.Data());
    return kFALSE;
  }
  char title[100];
  char name[100];
  if(hsBeta==0) {
    sprintf(title,"Cut: %.2f < #beta < %.2f;#beta",betaMin,betaMax);
    hsBeta= new TH1F("hsBeta",title,300,0.,1.5);
  }
  if(plStat==0) plStat= new TH2F("plStat",
    "Statistic;Sector of first tof;Sector of second tof",6,0.,6., 6,0.,6.);
  for(Int_t s=0;s<6;s++) for(Int_t m=0;m<4;m++) if(mdcsetup[s][m]!=0 && hsDCell[s][m]==0) {
    sprintf(name,"hsDCellS%iM%i",s+1,m+1);
    sprintf(title,"S.%i M.%i Cut: #Delta Cell <= %.1f;#Delta Cell",s+1,m+1,dCellWind[m]);
    hsDCell[s][m] = new TH1F(name,title,300,-30.,30.);
    for(Int_t l=0;l<6;l++) {
      sprintf(name,"hsDCellS%iM%iL%i",s+1,m+1,l+1);
      sprintf(title,"S.%i M.%i L.%i Cut: #Delta Cell <= %.1f;#Delta Cell",s+1,m+1,l+1,dCellWind[m]);
      hsDCellL[s][m][l] = new TH1F(name,title,300,-30.,30.);
    }
  }
  return kTRUE;
}

Bool_t HMdcCosmicsCandidate::finalize(void) {
  printf("\n-------------------------------------------------\n");
  if(nEventsSel == nEventsSelSc) {
    for(Int_t i=0;i<9;i++) if(nSelTr[i] > 0) {
      printf("%i modules: %10i tracks (%7.2f%c)\n",
      i,nSelTr[i],Float_t(nSelTr[i])/Float_t(nEventsSel)*100.,'%');
    }
    printf("  %i events filtered from %i.\n",nEventsSel,nEventsTot);
  } else {
    for(Int_t i=0;i<9;i++) if(nSelTr[i] > 0) {
      printf("%i modules: %10i tracks (%7.2f%c) ==> scaled down on %7.2f%c => %10i tracks (%7.2f%c)\n",
      i,nSelTr[i],Float_t(nSelTr[i])/Float_t(nEventsSel)*100.,'%',trackScaling[i]*100.,'%',
        nSelTrSc[i],Float_t(nSelTrSc[i])/Float_t(nEventsSelSc)*100.,'%');
    }
    printf("  %i events filtered from %i.\n",nEventsSelSc,nEventsTot);
  }
  printf("-------------------------------------------------\n");
  if(makeHists) {
    printf("  Histograms are stored to the file %s\n",histFileName.Data());
    // writing hists to the root file
    histFile->cd();
    if(hsBeta!=NULL) hsBeta->Write();
    if(plStat!=NULL) plStat->Write();
    for(Int_t s=0;s<6;s++) for(Int_t m=0;m<4;m++) if(hsDCell[s][m]!=NULL) {
      hsDCell[s][m]->Write();
      for(Int_t l=0;l<6;l++) hsDCellL[s][m][l]->Write();
    }
    histFile->Close();
  }
  return kTRUE;
}


void HMdcCosmicsCandidate::scaleDownTrack(UInt_t nMdcTr,Double_t sc) {
  if(nMdcTr<8) trackScaling[nMdcTr] = sc;
}

void HMdcCosmicsCandidate::scaleDownTrack(Double_t *sca) {
  for(Int_t nMdcs=0;nMdcs<9;nMdcs++) trackScaling[nMdcs] = sca[nMdcs];
}

HCosmicCalibEvSkip::HCosmicCalibEvSkip(void) 
    : HReconstructor("CosmicCalibEvSkip","CosmicCalibEvSkip") {
}

HCosmicCalibEvSkip::HCosmicCalibEvSkip(const Text_t *name,const Text_t *title) 
    : HReconstructor(name,title) {

}

Int_t HCosmicCalibEvSkip::execute(void) {
  return kSkipEvent;
}