/*
 * TofNucleiZCut.cpp
 *
 *  Created on: 21-apr-2009
 *      Author: Nicola Mori
 */

#ifndef NO_TOFNUCLEI

/*! @file TofNucleiZCut.cpp The TofNucleiZCut class implementation file */

#include "TofNucleiZCut.h"

const float TofNucleiZCut::_sigmaZdefault[] = { 0.12, 0.15, 0.20, 0.25, 0.32, 0.40, 0.45, 0.50 }; // the charge-width in a layer for Z=1..8
const unsigned int TofNucleiZCut::_maskArray[] = { S11, S12, S21, S22, S31, S32 }; //The layer - mask code association


TofNucleiZCut::TofNucleiZCut(const char *cutName, unsigned int Z, float lowerLimit, float upperLimit,
    unsigned int minLayers, unsigned int layersMask) :
  PamCut(cutName), _Z(Z), _mean(8, 6), _sigma(8, 6), _lowerLimit(lowerLimit), _upperLimit(upperLimit), _minLayers(
      minLayers), _layersMask(layersMask), _tofNuclei(NULL) {

  //TODO Put calibration parameters on a file; load only those parameters needed for the current Z.

  // Hard-coded calibrations
  // H
  _mean[0][0] = 0.987; // S11 peak
  _sigma[0][0] = 0.11; // S11 sigma
  _mean[0][1] = 0.988; // S12 peak
  _sigma[0][1] = 0.12; // S12 sigma
  _mean[0][2] = 0.961; // S21 peak
  _sigma[0][2] = 0.083;// S21 sigma
  _mean[0][3] = 0.966; // S22 peak
  _sigma[0][3] = 0.10; // S22 sigma
  _mean[0][4] = 0.981; // S31 peak
  _sigma[0][4] = 0.092;// S31 sigma
  _mean[0][5] = 0.979; // S32 peak
  _sigma[0][5] = 0.095;// S32 sigma

  // He
  _mean[1][0] = 1.96; // S11 peak
  _sigma[1][0] = 0.12;// S11 sigma
  _mean[1][1] = 1.95; // S12 peak
  _sigma[1][1] = 0.14;// S12 sigma
  _mean[1][2] = 1.95; // S21 peak
  _sigma[1][2] = 0.13;// S21 sigma
  _mean[1][3] = 1.96; // S22 peak
  _sigma[1][3] = 0.14;// S22 sigma
  _mean[1][4] = 1.99; // S31 peak
  _sigma[1][4] = 0.14;// S31 sigma
  _mean[1][5] = 2.00; // S32 peak
  _sigma[1][5] = 0.15;// S32 sigma

  // Heavier nuclei are not calibrated yet; using default values
  for (unsigned int z = 3; z < 9; z++) {
    for (int i = 0; i < 6; i++) {
      _mean[z - 1][i] = z;
      _sigma[z - 1][i] = _sigmaZdefault[z - 1];
    }
  }

#ifdef DEBUGPAMCUT

  TString hId;
  TString hTitle;

  for (UInt_t j = 0; j < 12; j++) {
    hId.Form("h_tof_z_vs_beta_%i", j);
    hTitle.Form("TOF Z vs beta (%i)", j);
    h_tof_z_beta[j] = new TH2F(hId.Data(), hTitle.Data(), 50, 0, 1, 50, 0, 10);
  }

#endif

}

int TofNucleiZCut::Check(PamLevel2 *event) {

  if (_Z < 1 || _Z > 8)
    return TOFNUCLEIZ_ILLEGALZ;

  // Check if ToFNuclei has already been initialized
  if (_tofNuclei == NULL)
    _tofNuclei = new ToFNuclei(event);

  Float_t *charge = _tofNuclei->Get_Charge_ToF_std_layer(); //Returns a 6 elements array
  unsigned int nGood = 0;
  bool good[6];

#ifdef DEBUGPAMCUT
  Int_t seqno = event->GetTrack(0)->GetTrkTrack()->GetSeqNo();
  Int_t ntrack_tof = 0;
  for (Int_t i=0; i<event->GetToFLevel2()->ntrk(); i++) {
    if ( event->GetToFLevel2()->GetToFTrkVar(i)->trkseqno == seqno ) ntrack_tof = i;
  }
  Float_t beta=event->GetToFLevel2()->CalcBeta(ntrack_tof,10.,10.,20.); // beta[12]
#endif

#ifdef DEBUGPAMCUT
  for (int i = 0; i < 6; i++) {
    h_tof_z_beta[i]->Fill(beta, charge[i]);
  }
#endif

  for (int i = 0; i < 6; i++) { // determine good, non-masked layers
    if ((charge[i] < 1000) && (charge[i] > 0) && ((_maskArray[i] & _layersMask) != _maskArray[i])) {
      good[i] = true;
      nGood++;
    }
    else {
      good[i] = false;
    }
  }

  if (nGood < _minLayers) {
    return TOFNUCLEIZ_TOOFEWLAYERS;
  }

  for (int i = 0; i < 6; i++) { // check done only for good, non-masked layers
    if (good[i]) {
      if ((charge[i] > _mean[_Z - 1][i] + _upperLimit * _sigma[_Z - 1][i]) || ((charge[i] < _mean[_Z - 1][i]
          - _lowerLimit * _sigma[_Z - 1][i]))) {
        return TOFNUCLEIZ_OUTOFBOUNDS;
      }
    }
  }

#ifdef DEBUGPAMCUT
  for (int i = 0; i < 6; i++) {
    h_tof_z_beta[6+i]->Fill(beta, charge[i]);
  }
#endif

  return CUTOK;

}

#endif /* NO_TOFNUCLEI */