ToFLevel2/src/ToFLevel2.cpp

/**
 * \file ToFLevel2.cpp
 * \author Gianfranca DeRosa, Wolfgang Menn
 *
 * WM dec 2008: Description of "GetdEdx" changed
 * WM dec 2008: "GetdEdxPaddle" modified: Now includes saturation limit
 *              PMTs higher than the saturation limit are not used for dEdx
 */

#include <ToFLevel2.h>
using namespace std;
ClassImp(ToFPMT);
ClassImp(ToFTrkVar);
ClassImp(ToFLevel2);

ToFPMT::ToFPMT(){
  pmt_id = 0;
  adc = 0.;
  tdc_tw = 0.;
  tdc = 0.;
}

ToFPMT::ToFPMT(const ToFPMT &t){
  pmt_id = t.pmt_id;
  adc = t.adc;
  tdc_tw = t.tdc_tw;
  tdc = t.tdc;
}

void ToFPMT::Clear(Option_t *t){
  pmt_id = 0;
  adc = 0.;
  tdc_tw = 0.;
  tdc = 0.;
}


ToFTrkVar::ToFTrkVar() {
  trkseqno = 0;
  npmttdc = 0;
  npmtadc = 0;
  pmttdc = TArrayI(48);
  pmtadc = TArrayI(48);
  tdcflag = TArrayI(48); // gf: 30 Nov 2006
  adcflag = TArrayI(48); // gf: 30 Nov 2006
  dedx = TArrayF(48);
  //
  //
  memset(beta,  0, 13*sizeof(Float_t));
  memset(xtofpos,  0, 3*sizeof(Float_t));
  memset(ytofpos,  0, 3*sizeof(Float_t));
  memset(xtr_tof,  0, 6*sizeof(Float_t));
  memset(ytr_tof,  0, 6*sizeof(Float_t));
  //
};

void ToFTrkVar::Clear(Option_t *t) {
  trkseqno = 0;
  npmttdc = 0;
  npmtadc = 0;
  pmttdc.Reset();
  pmtadc.Reset();
  tdcflag.Reset(); // gf: 30 Nov 2006
  adcflag.Reset(); // gf: 30 Nov 2006
  dedx.Reset(); 
  //
  memset(beta,  0, 13*sizeof(Float_t));
  memset(xtofpos,  0, 3*sizeof(Float_t));
  memset(ytofpos,  0, 3*sizeof(Float_t));
  memset(xtr_tof,  0, 6*sizeof(Float_t));
  memset(ytr_tof,  0, 6*sizeof(Float_t));
  //
};

ToFTrkVar::ToFTrkVar(const ToFTrkVar &t){

  trkseqno = t.trkseqno;  
  //
  npmttdc = t.npmttdc;
  npmtadc = t.npmtadc;
  (t.pmttdc).Copy(pmttdc);
  (t.pmtadc).Copy(pmtadc);
  (t.tdcflag).Copy(tdcflag); // gf: 30 Nov 2006
  (t.adcflag).Copy(adcflag); // gf: 30 Nov 2006
  (t.dedx).Copy(dedx);
  //
  memcpy(beta,t.beta,sizeof(beta));
  memcpy(xtofpos,t.xtofpos,sizeof(xtofpos));
  memcpy(ytofpos,t.ytofpos,sizeof(ytofpos));
  memcpy(xtr_tof,t.xtr_tof,sizeof(xtr_tof));
  memcpy(ytr_tof,t.ytr_tof,sizeof(ytr_tof));
  //
};

ToFLevel2::ToFLevel2() {    
  //
//  PMT = new TClonesArray("ToFPMT",12); //ELENA
//  ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
  PMT = 0; //ELENA
  ToFTrk = 0; //ELENA
  //
  this->Clear();
  //
};

void ToFLevel2::Set(){//ELENA
    if(!PMT)PMT = new TClonesArray("ToFPMT",12); //ELENA
    if(!ToFTrk)ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
}//ELENA

void ToFLevel2::Clear(Option_t *t){
  //
  if(ToFTrk)ToFTrk->Delete(); //ELENA
  if(PMT)PMT->Delete(); //ELENA
  memset(tof_j_flag, 0, 6*sizeof(Int_t));
  unpackError = 0;
  // 
};

void ToFLevel2::Delete(Option_t *t){ //ELENA
  //
  if(ToFTrk){ 
      ToFTrk->Delete(); //ELENA
      delete ToFTrk;  //ELENA
  }
  if(PMT){
      PMT->Delete(); //ELENA
      delete PMT; //ELENA
  } //ELENA
  // 
}; //ELENA

ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t itrk){
  //    
  if(itrk >= ntrk()){
    printf(" ToFLevel2 ERROR: track related variables set %i does not exists! \n",itrk);
    printf("                      stored track related variables = %i \n",ntrk());
    return(NULL);
  }  
  //
  if(!ToFTrk)return 0; //ELENA
  TClonesArray &t = *(ToFTrk);
  ToFTrkVar *toftrack = (ToFTrkVar*)t[itrk];
  return toftrack;
}

ToFPMT *ToFLevel2::GetToFPMT(Int_t ihit){
  //    
  if(ihit >= npmt()){
    printf(" ToFLevel2 ERROR: pmt variables set %i does not exists! \n",ihit);
    printf("                  stored pmt variables = %i \n",npmt());
    return(NULL);
  }  
  //
  if(!PMT)return 0; //ELENA
  TClonesArray &t = *(PMT);
  ToFPMT *tofpmt = (ToFPMT*)t[ihit];
  return tofpmt;
}
//--------------------------------------
//
// 
//--------------------------------------
/**
 * Method to get the plane ID (11 12 21 22 31 32) from the plane index (0 1 2 3 4 5)
 * @param Plane index (0,1,2,3,4,5).
 */
  Int_t  ToFLevel2::GetToFPlaneID(Int_t ip){ 
      if(ip>=0 && ip<6)return 10*((int)(ip/2+1.1))+(ip%2)+1;
      else return -1;
  };
/**
 * Method to get the plane index (0 1 2 3 4 5) from the plane ID (11 12 21 22 31 32)
 * @param plane Plane ID (11, 12, 21, 22, 31, 32) 
 */
  Int_t  ToFLevel2::GetToFPlaneIndex(Int_t plane_id){ 
      if(
         plane_id == 11 ||
         plane_id == 12 ||
         plane_id == 21 ||
         plane_id == 22 ||
         plane_id == 31 ||
         plane_id == 32 ||
         false)return (Int_t)(plane_id/10)*2-1- plane_id%2;
      else return -1;
  };
/**
 * Method to know if a given ToF paddle was hit, that is there is a TDC signal 
 * from both PMTs. The method uses the "tof_j_flag" variable.
 * @param plane Plane ID (11, 12, 21, 22, 31, 32) or Plane index (0,1,2,3,4,5).
 * @param paddle_id Paddle ID.
 * @return 1 if the paddle was hit.
 */
  Bool_t ToFLevel2::HitPaddle(Int_t plane, Int_t paddle_id){  //<<< NEW
    Int_t ip = -1;
    if     (plane>=6             ) ip = GetToFPlaneIndex(plane);
    else if(plane>=0 && plane < 6) ip = plane;
    Int_t flag=0;
    if(ip != -1)flag = tof_j_flag[ip] & (int)pow(2.,(double)paddle_id);
    if(
       (ip == 0 && paddle_id < 8 && flag) ||
       (ip == 1 && paddle_id < 6 && flag) ||
       (ip == 2 && paddle_id < 2 && flag) ||
       (ip == 3 && paddle_id < 2 && flag) ||
       (ip == 4 && paddle_id < 3 && flag) ||
       (ip == 5 && paddle_id < 3 && flag) ||
       false) return true;
    else return false;
};
/**
 * Method to get the number of hit paddles on a ToF plane.
 * @param plane Plane ID (11, 12, 21, 22, 31, 32) or Plane index (0,1,2,3,4,5).
 */
Int_t ToFLevel2::GetNHitPaddles(Int_t plane){
    Int_t npad=0;
    for(Int_t i=0; i<8; i++)npad = npad + (int)HitPaddle(plane,i);
    return npad;
};

//wm Nov 08
//gf Apr 07
/**
 * Method to get the mean dEdx from a ToF layer - ATTENTION:
 * It will sum up the dEdx of all the paddles, but since by definition
 * only the paddle hitted by the track gets a dEdx value and the other 
 * paddles are set to zero, the output is just the dEdx of the hitted 
 * paddle in each layer! 
 * The "adcfl" option is not very useful (an artificial dEdx is per
 * definition= 1 mip and not a real measurement), anyway left in the code 
 * @param notrack Track Number
 * @param plane Plane index (0,1,2,3,4,5)
 * @param adcflag in the plane (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )
 */
Float_t ToFLevel2::GetdEdx(Int_t notrack, Int_t plane, Int_t adcfl){

  Float_t dedx = 0.;
  Float_t PadEdx =0.;
  Int_t SatWarning;
  Int_t pad=-1;
  //
  ToFTrkVar *trk = GetToFTrkVar(notrack);
  if(!trk) return 0; //ELENA
  //
  for (Int_t ii=0; ii<GetNPaddle(plane); ii++){
    Int_t paddleid=ii;
    pad = GetPaddleid(plane,paddleid);
    GetdEdxPaddle(notrack, pad, adcfl, PadEdx, SatWarning);
    dedx += PadEdx; 
  };
  //
  return(dedx);
};

/**
 * Method to fill the ADC_C 4x12 matrix with the dEdx values and the TDC 4x12 matrix 
 * with the time-walk corrected TDC values.
 * @param notrack Track Number
 * @param adc  ADC_C matrix with dEdx values
 * @param tdc  TDC matrix
 */
void ToFLevel2::GetMatrix(Int_t notrack, Float_t adc[4][12], Float_t tdc[4][12]){
  //
  for (Int_t aa=0; aa<4;aa++){
    for (Int_t bb=0; bb<12;bb++){
      adc[aa][bb] = 1000.;
      tdc[aa][bb] = 4095.;
    };
  };
  //
  Int_t pmt_id = 0;
  Int_t hh = 0;
  Int_t kk = 0;
  //
  ToFTrkVar *trk = GetToFTrkVar(notrack);
  if(!trk)return; //ELENA
  //
  for (Int_t i=0; i<trk->npmtadc; i++){
    //
    pmt_id = (trk->pmtadc).At(i);
    //
    GetPMTIndex(pmt_id,hh,kk);
    adc[kk][hh] = (trk->dedx).At(i);   
    //
  };
  //
  for (Int_t i=0; i<npmt(); i++){
    //
    ToFPMT *pmt = GetToFPMT(i);
    if(!pmt)break; //ELENA
    //
    GetPMTIndex(pmt->pmt_id,hh,kk);
    //
    tdc[kk][hh] = pmt->tdc_tw;   
    //
  };
  //
  return;
};


/**
 * Method to get the plane index (0 - 5) for the PMT_ID as input
 * @param pmt_id  PMT_ID (0 - 47)
 */
Int_t ToFLevel2::GetPlaneIndex(Int_t pmt_id){
  TString pmtname = GetPMTName(pmt_id);
  pmtname.Resize(3);
  if ( !strcmp(pmtname,"S11") ) return(0);
  if ( !strcmp(pmtname,"S12") ) return(1);
  if ( !strcmp(pmtname,"S21") ) return(2);
  if ( !strcmp(pmtname,"S22") ) return(3);
  if ( !strcmp(pmtname,"S31") ) return(4);
  if ( !strcmp(pmtname,"S32") ) return(5);
  return(-1);
};


/**
 * Method to get the PMT_ID if the index (4,12) is given. We have 4 channels on
 * each of the 12 half-boards, this method decodes which PMT is cables to which 
 * channel.
 * @param hh Channel
 * @param kk HalfBoard
 */
Int_t ToFLevel2::GetPMTid(Int_t hh, Int_t kk){
  //
  short tof[4][24] = {
    {4, 4,  4,  4,  1,  1, 2, 2,  3,  3, 3, 3,  3,  3, 1, 1,  1,  1, 2, 3,  3, 3, 3,  4},
    {1, 3,  5,  7, 10, 12, 2, 4,  2,  4, 6, 8, 10, 12, 1, 5,  3,  9, 7, 9, 11, 1, 5,  9},
    {2, 2,  2,  2,  1,  1, 1, 1,  4,  4, 4, 4,  4,  4, 2, 1,  2,  1, 2, 2,  2, 3, 3,  4},
    {6, 8, 12, 10,  8,  6, 4, 2, 12, 10, 8, 6,  4,  2, 9, 7, 11, 11, 5, 3,  1, 3, 7, 11}
  };
  //
  Int_t ind = 0;
  Int_t k = 0;
  while (k < 24){
    Int_t j = 0;
    while (j < 2){
      Int_t ch = tof[2*j][k]     - 1;
      Int_t hb = tof[2*j + 1][k] - 1;      
      /* tofEvent->tdc[ch][hb] */      
      if( ch == hh && hb == kk ){       
        ind = 2*k + j;
        break;
      };
      j++;
    };
    k++;
  };
  return ind;
};


/**
 * Method to get the PMT index if the PMT ID is given. This method is the
 * "reverse" of method "GetPMTid"
 * @param ind  PMT_ID (0 - 47)
 * @param hb   HalfBoard
 * @param ch   Channel
 */
void ToFLevel2::GetPMTIndex(Int_t ind, Int_t &hb, Int_t &ch){
  //
  short tof[4][24] = {
    {4, 4,  4,  4,  1,  1, 2, 2,  3,  3, 3, 3,  3,  3, 1, 1,  1,  1, 2, 3,  3, 3, 3,  4},
    {1, 3,  5,  7, 10, 12, 2, 4,  2,  4, 6, 8, 10, 12, 1, 5,  3,  9, 7, 9, 11, 1, 5,  9},
    {2, 2,  2,  2,  1,  1, 1, 1,  4,  4, 4, 4,  4,  4, 2, 1,  2,  1, 2, 2,  2, 3, 3,  4},
    {6, 8, 12, 10,  8,  6, 4, 2, 12, 10, 8, 6,  4,  2, 9, 7, 11, 11, 5, 3,  1, 3, 7, 11}
  };
  //
  Int_t k = 0;
  while (k < 24){
    Int_t j = 0;
    while (j < 2){
      /* tofEvent->tdc[ch][hb] */            
      if( ind == 2*k + j ){
        ch = tof[2*j][k]     - 1;
        hb = tof[2*j + 1][k] - 1;      
        return;
      };
      j++;
    };
    k++;
  };
  return;
};


//  wm Nov 08 revision - saturation values included
/// gf Apr 07
/**
 * Method to get the dEdx from a given ToF paddle.
 * If two PMTs are good, the mean dEdx of both PMTs is taken, otherwise
 * just the dEdx of the "good" PMT. If both PMTs are above saturation => dEdx=1000 
 * @param notrack Track Number
 * @param Paddle index (0,1,...,23).
 * @param adcflag in the paddle (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )
 * @param PadEdx dEdx from a given ToF paddle
 * @param SatWarning 1 if the PMT ios near saturation region (adcraw ~3000)
 */
void ToFLevel2::GetdEdxPaddle(Int_t notrack, Int_t paddleid, Int_t adcfl, Float_t &PadEdx, Int_t &SatWarning){

/*
Float_t  PMTsat[48] = {
3162.14, 3165.48, 3153.85, 3085.73, 3089.65, 3107.64, 3097.52, 3078.37, 
3130.05, 3087.07, 3112.22, 3102.92, 3080.58, 3092.55, 3087.94, 3125.03, 
3094.09, 3143.16, 3125.51, 3181.27, 3092.09, 3124.98, 3069.3, 3095.53, 
3097.11, 3133.53, 3114.73, 3113.01, 3091.19, 3097.99, 3033.84, 3134.98, 
3081.37, 3111.04, 3066.77, 3108.17, 3133, 3111.06, 3052.52, 3140.66, 
3106.33, 3094.85, 3150.85, 3118.8, 3096.24, 3118.47,3111.36, 3117.11 } ;
*/

// new values from Napoli dec 2008
Float_t  PMTsat[48] = {
3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };

for (Int_t i=0; i<48;i++) PMTsat[i] = PMTsat[i] - 5.;  // safety margin 


  PadEdx = 0.;
//  SatWarning = 1000;
  SatWarning = 0;   // 0=good, increase for each bad PMT

  Float_t dEdx[48] = {0};
  Int_t pmt_id = -1;
  Float_t adcraw[48];
  //
  ToFTrkVar *trk = GetToFTrkVar(notrack);
  if(!trk) return; //ELENA
  //

  Int_t pmtleft=-1;
  Int_t pmtright=-1;
  GetPaddlePMT(paddleid, pmtleft, pmtright);

  adcraw[pmtleft] = 4095;
  adcraw[pmtright] = 4095;

  
  for (Int_t jj=0; jj<npmt(); jj++){
    
    ToFPMT *pmt = GetToFPMT(jj);
    if(!pmt)break; //ELENA
    
    pmt_id = pmt->pmt_id;
    if(pmt_id==pmtleft){
      adcraw[pmtleft] = pmt->adc;
    }
    
    if(pmt_id==pmtright){
      adcraw[pmtright] = pmt->adc;
    }
  }

  
  for (Int_t i=0; i<trk->npmtadc; i++){

    if((trk->adcflag).At(i)==0 || adcfl==100){
      if((trk->pmtadc).At(i) == pmtleft)dEdx[pmtleft] = (trk->dedx).At(i);
      if((trk->pmtadc).At(i) == pmtright)dEdx[pmtright] = (trk->dedx).At(i);
    }else{ 
      if((trk->pmtadc).At(i) == pmtleft)dEdx[pmtleft] = 0.;
      if((trk->pmtadc).At(i) == pmtright)dEdx[pmtright] = 0.;
    }
  }


//  if( adcraw[pmtleft] >3000 || adcraw[pmtright] >3000)SatWarning=1;  //old version

// Increase SatWarning Counter for each PMT>Sat
  if( adcraw[pmtleft] > PMTsat[pmtleft])SatWarning++;  
  if( adcraw[pmtright] > PMTsat[pmtright])SatWarning++;

// if ADC  > sat set dEdx=1000
  if( adcraw[pmtleft] > PMTsat[pmtleft]) dEdx[pmtleft] = 1000.;
  if( adcraw[pmtright] > PMTsat[pmtright]) dEdx[pmtright] = 1000. ;

// if two PMT are good, take mean dEdx, otherwise only the good dEdx
  if(dEdx[pmtleft]<1000 && dEdx[pmtright]<1000) PadEdx = (dEdx[pmtleft]+dEdx[pmtright])*0.5;
  if(dEdx[pmtleft]==1000 && dEdx[pmtright]<1000) PadEdx = dEdx[pmtright];  
  if(dEdx[pmtleft]<1000 && dEdx[pmtright]==1000) PadEdx = dEdx[pmtleft];
  
};
//


// gf Apr 07

/**
 * Method to get the PMT name (like "S11_1A") if the PMT_ID is given.
 * Indexes of corresponding  plane, paddle and  pmt are also given as output.
 * @param ind  PMT_ID (0 - 47)
 * @param iplane plane index (0 - 5)
 * @param ipaddle paddle index (relative to the plane)
 * @param ipmt pmt index (0(A), 1(B))
 */
TString ToFLevel2::GetPMTName(Int_t ind, Int_t &iplane, Int_t &ipaddle,Int_t &ipmt){
  
  TString pmtname = " ";
  
  TString photoS[48] = {
    "S11_1A", "S11_1B", "S11_2A", "S11_2B", "S11_3A", "S11_3B", "S11_4A",
    "S11_4B",
    "S11_5A", "S11_5B", "S11_6A", "S11_6B", "S11_7A", "S11_7B", "S11_8A",
    "S11_8B",
    "S12_1A", "S12_1B", "S12_2A", "S12_2B", "S12_3A", "S12_3B", "S12_4A",
    "S12_4B", "S12_5A",  "S12_5B", "S12_6A", "S12_6B",
    "S21_1A", "S21_1B", "S21_2A", "S21_2B",
    "S22_1A", "S22_1B", "S22_2A", "S22_2B",
    "S31_1A", "S31_1B", "S31_2A", "S31_2B", "S31_3A", "S31_3B",
    "S32_1A", "S32_1B", "S32_2A", "S32_2B", "S32_3A", "S32_3B"
  };
  
  
  pmtname = photoS[ind].Data();
  
  TString ss = pmtname(1,2);
  iplane  = (int)(atoi(ss.Data())/10)*2-3+atoi(ss.Data())%10;
  ss = pmtname(4);
  ipaddle = atoi(ss.Data())-1 ;
  if( pmtname.Contains("A") )ipmt=0;
  if( pmtname.Contains("B") )ipmt=1;
  
  return pmtname;
};
/**
 * Method to get the PMT name (like "S11_1A") if the PMT_ID is given
 * @param ind  PMT_ID (0 - 47)
 */
TString ToFLevel2::GetPMTName(Int_t ind){

  Int_t iplane  = -1;
  Int_t ipaddle = -1;
  Int_t ipmt    = -1;
  return GetPMTName(ind,iplane,ipaddle,ipmt);
  
};

// wm jun 08
Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane){
return GetPaddleIdOfTrack(xtr ,ytr ,plane, 0.4);
}

// gf Apr 07
Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane, Float_t margin){
  
  Double_t xt,yt,xl,xh,yl,yh;
  
  Float_t tof11_x[8] = {-17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85};
  Float_t tof12_y[6] = { -13.75,-8.25,-2.75,2.75,8.25,13.75};
  Float_t tof21_y[2] = { 3.75,-3.75};
  Float_t tof22_x[2] = { -4.5,4.5};
  Float_t tof31_x[3] = { -6.0,0.,6.0};
  Float_t tof32_y[3] = { -5.0,0.0,5.0};
  
  //  S11 8 paddles  33.0 x 5.1 cm
  //  S12 6 paddles  40.8 x 5.5 cm
  //  S21 2 paddles  18.0 x 7.5 cm
  //  S22 2 paddles  15.0 x 9.0 cm
  //  S31 3 paddles  15.0 x 6.0 cm
  //  S32 3 paddles  18.0 x 5.0 cm
  
  Int_t paddleidoftrack=-1;
  //
  
  //--- S11 ------
  
  if(plane==0){
    xt = xtr;
    yt = ytr;
    paddleidoftrack=-1;
    yl = -33.0/2. ;
    yh =  33.0/2. ;
    if ((yt>yl)&&(yt<yh)) {
      for (Int_t i1=0; i1<8;i1++){
        xl = tof11_x[i1] - (5.1-margin)/2. ;
        xh = tof11_x[i1] + (5.1-margin)/2. ;
        if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
      }
    }
  }
  //      cout<<"S11  "<<paddleidoftrack[0]<<"\n";
  
  //--- S12 -------
  if(plane==1){
    xt = xtr;
    yt = ytr;
    paddleidoftrack=-1;
    xl = -40.8/2. ;
    xh =  40.8/2. ;
    
    if ((xt>xl)&&(xt<xh)) {
      for (Int_t i1=0; i1<6;i1++){
        yl = tof12_y[i1] - (5.5-margin)/2. ;
        yh = tof12_y[i1] + (5.5-margin)/2. ;
        if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;
      }
    }
  }
  
  //--- S21 ------

  if(plane==2){
    xt = xtr;
    yt = ytr;
    paddleidoftrack=-1;
    xl = -18./2. ;
    xh =  18./2. ;
    
    if ((xt>xl)&&(xt<xh)) {
      for (Int_t i1=0; i1<2;i1++){
        yl = tof21_y[i1] - (7.5-margin)/2. ;
        yh = tof21_y[i1] + (7.5-margin)/2. ;
        if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;
      }
    }
  }
  
  //--- S22 ------
  if(plane==3){
    xt = xtr;
    yt = ytr;
    paddleidoftrack=-1;
    yl = -15./2. ;
    yh =  15./2. ;
    
    if ((yt>yl)&&(yt<yh)) {
      for (Int_t i1=0; i1<2;i1++){
        xl = tof22_x[i1] - (9.0-margin)/2. ;
        xh = tof22_x[i1] + (9.0-margin)/2. ;
        if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
      }
    }
  }  
  
  //--- S31 ------
  if(plane==4){
    xt = xtr;
    yt = ytr;
    paddleidoftrack=-1;
    yl = -15.0/2. ;
    yh =  15.0/2. ;
    
    if ((yt>yl)&&(yt<yh)) {
      for (Int_t i1=0; i1<3;i1++){
        xl = tof31_x[i1] - (6.0-margin)/2. ;
        xh = tof31_x[i1] + (6.0-margin)/2. ;
        if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
      }
    }
  }  
  
  //---  S32 ------
  if(plane==5){
    xt = xtr;
    yt = ytr;
    paddleidoftrack=-1;
    xl = -18.0/2. ;
    xh =  18.0/2. ;
    
    if ((xt>xl)&&(xt<xh)) {
      for (Int_t i1=0; i1<3;i1++){
        yl = tof32_y[i1] - (5.0-margin)/2. ;
        yh = tof32_y[i1] + (5.0-margin)/2. ;
        if ((yt>yl)&&(yt<yh)) paddleidoftrack=i1;
      }
    }
  }
  
  return paddleidoftrack;

}  

//

// gf Apr 07

void ToFLevel2::GetPMTPaddle(Int_t pmt_id, Int_t &plane, Int_t &paddle){
  
  plane = GetPlaneIndex(pmt_id);

  if(plane == 0){
    if(pmt_id==0 || pmt_id==1)paddle=0;
    if(pmt_id==2 || pmt_id==3)paddle=1;
    if(pmt_id==4 || pmt_id==5)paddle=2;
    if(pmt_id==6 || pmt_id==7)paddle=3;
    if(pmt_id==8 || pmt_id==9)paddle=4;
    if(pmt_id==10 || pmt_id==11)paddle=5;
    if(pmt_id==12 || pmt_id==13)paddle=6;
    if(pmt_id==14 || pmt_id==15)paddle=7;
  }
  
  if(plane == 1){
    if(pmt_id==16 || pmt_id==17)paddle=0;
    if(pmt_id==18 || pmt_id==19)paddle=1;
    if(pmt_id==20 || pmt_id==21)paddle=2;
    if(pmt_id==22 || pmt_id==23)paddle=3;
    if(pmt_id==24 || pmt_id==25)paddle=4;
    if(pmt_id==26 || pmt_id==27)paddle=5;
  }
  
  if(plane == 2){
    if(pmt_id==28 || pmt_id==29)paddle=0;
    if(pmt_id==30 || pmt_id==31)paddle=1;
  }
  
  if(plane == 3){
    if(pmt_id==32 || pmt_id==33)paddle=0;
    if(pmt_id==34 || pmt_id==35)paddle=1;
  }
  
  if(plane == 4){
    if(pmt_id==36 || pmt_id==37)paddle=0;
    if(pmt_id==38 || pmt_id==39)paddle=1;
    if(pmt_id==40 || pmt_id==41)paddle=2;
  }
  
  if(plane == 5){
    if(pmt_id==42 || pmt_id==43)paddle=0;
    if(pmt_id==44 || pmt_id==45)paddle=1;
    if(pmt_id==46 || pmt_id==47)paddle=2;
  }
  return;
}

//

// gf Apr 07

void ToFLevel2::GetPaddlePMT(Int_t paddle, Int_t &pmtleft, Int_t &pmtright){
  pmtleft=paddle*2;
  pmtright= pmtleft+1;  
  return;
}

//


// // gf Apr 07

void ToFLevel2::GetPaddleGeometry(Int_t plane, Int_t paddle, Float_t &xleft, Float_t &xright, Float_t &yleft, Float_t &yright){
  
  Int_t i1;

  Float_t tof11_x[8] = {-17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85};
  Float_t tof12_y[6] = { -13.75,-8.25,-2.75,2.75,8.25,13.75};
  Float_t tof21_y[2] = { 3.75,-3.75};
  Float_t tof22_x[2] = { -4.5,4.5};
  Float_t tof31_x[3] = { -6.0,0.,6.0};
  Float_t tof32_y[3] = { -5.0,0.0,5.0};
        
  //  S11 8 paddles  33.0 x 5.1 cm
  //  S12 6 paddles  40.8 x 5.5 cm
  //  S21 2 paddles  18.0 x 7.5 cm
  //  S22 2 paddles  15.0 x 9.0 cm
  //  S31 3 paddles  15.0 x 6.0 cm
  //  S32 3 paddles  18.0 x 5.0 cm

  if(plane==0)
    {
      for (i1=0; i1<8;i1++){
        if(i1 == paddle){
          xleft = tof11_x[i1] - 5.1/2.;
          xright = tof11_x[i1] + 5.1/2.;
          yleft = -33.0/2.;
          yright = 33.0/2.;
        }
      }
    }
  
  if(plane==1)
    {
      for (i1=0; i1<6;i1++){
        if(i1 == paddle){
          xleft = -40.8/2.;
          xright = 40.8/2.;
          yleft = tof12_y[i1] - 5.5/2.;
          yright = tof12_y[i1] + 5.5/2.;
        }
      }
    }

  if(plane==2)
    {
      for (i1=0; i1<2;i1++){
        if(i1 == paddle){
          xleft =  -18./2.;
          xright = 18./2.;
          yleft = tof21_y[i1] - 7.5/2.;
          yright = tof21_y[i1] + 7.5/2.;
        }
      }
    }
  
  if(plane==3)
    {
      for (i1=0; i1<2;i1++){
        if(i1 == paddle){
          xleft = tof22_x[i1] - 9.0/2.;
          xright = tof22_x[i1] + 9.0/2.;
          yleft = -15./2.;
          yright = 15./2.;
        }
      }
    }


  if(plane==4)
    {
      for (i1=0; i1<3;i1++){
        if(i1 == paddle){
          xleft = tof31_x[i1] - 6.0/2.;
          xright = tof31_x[i1] + 6.0/2.;
          yleft = -15./2.;
          yright = 15./2.;
        }
      }
    }

  if(plane==5)
    {
      for (i1=0; i1<3;i1++){
        if(i1 == paddle){
          xleft = -18.0/2.;
          xright = 18.0/2.;
          yleft = tof32_y[i1] - 5.0/2.;
          yright = tof32_y[i1] + 5.0/2.;
        }
      }
    }
  return;
}

// gf Apr 07
/**
 * Method to get the paddle index (0,...23) if the plane ID and the paddle id in the plane is given. 
 * This method is the
 * "reverse" of method "GetPaddlePlane"
 * @param plane    (0 - 5)
 * @param paddle   (plane=0, paddle = 0,...5)
 * @param padid    (0 - 23)
 */
Int_t ToFLevel2::GetPaddleid(Int_t plane, Int_t paddle)
{
  Int_t padid=-1;
  Int_t pads[6]={8,6,2,2,3,3};

  int somma=0;
  int np=plane;
  for(Int_t j=0; j<np; j++){
    somma+=pads[j];
  }
  padid=paddle+somma;
  return padid;

}


// gf Apr 07
/**
 * Method to get the plane ID and the paddle id in the plane if the paddle index (0,...23) is given. 
 * This method is the
 * "reverse" of method "GetPaddleid"
 * @param pad      (0 - 23)
 * @param plane    (0 - 5)
 * @param paddle   (plane=0, paddle = 0,...5)
 */
void ToFLevel2::GetPaddlePlane(Int_t pad, Int_t &plane, Int_t &paddle)
{

  Int_t pads11=8;
  Int_t pads12=6;
  Int_t pads21=2;
  Int_t pads22=2;
  Int_t pads31=3;
  // Int_t pads32=3;

  if(pad<8){
    plane=0;
    paddle=pad;
    return;
  }

  if(7<pad<14){
    plane=1;
    paddle=pad-pads11;
    return;
  }
  
  if(13<pad<16){
    plane=2;
    paddle=pad-pads11-pads12;
    return;
  }

  if(15<pad<18){
    plane=3;
    paddle=pad-pads11-pads12-pads21;
    return;
  }

  if(17<pad<21){
    plane=4;
    paddle=pad-pads11-pads12-pads21-pads22;
    return;
  }

  if(20<pad<24){
    plane=5;
    paddle=pad-pads11-pads12-pads21-pads22-pads31;
    return;
  }   

}


Int_t ToFLevel2::GetNPaddle(Int_t plane){

  Int_t npaddle=-1;

  Int_t pads11=8;
  Int_t pads12=6;
  Int_t pads21=2;
  Int_t pads22=2;
  Int_t pads31=3;
  Int_t pads32=3;

  if(plane==0)npaddle=pads11;
  if(plane==1)npaddle=pads12;
  if(plane==2)npaddle=pads21;
  if(plane==3)npaddle=pads22;
  if(plane==4)npaddle=pads31;
  if(plane==5)npaddle=pads32;

  return npaddle;

}


/// wm feb 08

/**
 * Method to calculate Beta from the 12 single measurements
 * we check the individual weights for artificial TDC values, then calculate
 * am mean beta for the first time. In a second step we loop again through
 * the single measurements, checking for the residual from the mean
 * The cut on the residual reject measurements > "x"-sigma. A chi2 value is
 * calculated, furthermore a "quality" value by adding the weights which
 * are finally used. If all measurements are taken, "quality" will be = 22.47.
 * A chi2 cut around 3-4 and a quality-cut > 20 is needed for clean beta
 * measurements like antiprotons etc.
 * The Level2 output is derived in the fortran routines using: 10.,10.,20.
 * @param notrack Track Number
 * @param cut on residual: difference between single measurement and mean 
 * @param cut on "quality"
 * @param cut on chi2
 */

Float_t ToFLevel2::CalcBeta(Int_t notrack, Float_t resmax, Float_t qualitycut, Float_t chi2cut){

//  cout<<" in CalcBeta "<<resmax<<" "<<chi2cut<<" "<<qualitycut<<endl;

  Float_t bxx = 100.;
  //
  ToFTrkVar *trk = GetToFTrkVar(notrack);
  if(!trk) return 0; //ELENA


  Float_t chi2,xhelp,beta_mean;
  Float_t w_i[12],quality,sw,sxw,res,betachi,beta_mean_inv;
  Float_t b[12],tdcfl;
  Int_t  pmt_id,pmt_plane;

  for (Int_t i=0; i<12; i++){
    b[i] = trk->beta[i];
                              }
      

//========================================================================
//---  Find out ToF layers with artificial TDC values & fill vector    ---
//========================================================================

Float_t  w_il[6];

     for (Int_t jj=0; jj<6;jj++) {
         w_il[jj] = 1000.;
                                 }


  for (Int_t i=0; i<trk->npmttdc; i++){
    //
    pmt_id = (trk->pmttdc).At(i);
    pmt_plane = GetPlaneIndex(pmt_id);
    tdcfl = (trk->tdcflag).At(i);
    if (w_il[pmt_plane] != 1.) w_il[pmt_plane] = tdcfl; //tdcflag
                                     };
  
//========================================================================
//---  Set weights for the 12 measurements using information for top and bottom:
//---  if no measurements: weight = set to very high value=> not used
//---  top or bottom artificial: weight*sqrt(2)
//---  top and bottom artificial: weight*sqrt(2)*sqrt(2)
//========================================================================

Int_t itop[12] = {0,0,1,1,2,2,3,3,0,0,1,1};
Int_t ibot[12] = {4,5,4,5,4,5,4,5,2,3,2,3};

     xhelp= 1E09;
   
     for (Int_t jj=0; jj<12;jj++) {
     if (jj<4)           xhelp = 0.11;    // S1-S3
     if ((jj>3)&&(jj<8)) xhelp = 0.18;    // S2-S3
     if (jj>7)           xhelp = 0.28;    // S1-S2
     if ((w_il[itop[jj]] == 1000.) && (w_il[ibot[jj]] == 1000.)) xhelp = 1E09;
     if ((w_il[itop[jj]] == 1) || (w_il[ibot[jj]] == 1.)) xhelp = xhelp*1.414 ;
     if ((w_il[itop[jj]] == 1) && (w_il[ibot[jj]] == 1.)) xhelp = xhelp*2. ;

     w_i[jj] = 1./xhelp;
                                  }


//========================================================================
//--- Calculate mean beta for the first time -----------------------------
//--- We are using "1/beta" since its error is gaussian ------------------
//========================================================================

      Int_t icount=0;
      sw=0.;
      sxw=0.;
      beta_mean=100.;

          for (Int_t jj=0; jj<12;jj++){
        if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.))
         {
            icount= icount+1;
            sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
            sw =sw + w_i[jj]*w_i[jj] ;

         }
         }

      if (icount>0) beta_mean=1./(sxw/sw);
      beta_mean_inv = 1./beta_mean;

//========================================================================
//--- Calculate beta for the second time, use residuals of the single
//--- measurements to get a chi2 value
//========================================================================

      icount=0;
      sw=0.;
      sxw=0.;
      betachi = 100.;
      chi2 = 0.;
      quality=0.;


          for (Int_t jj=0; jj<12;jj++){
       if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.)&&(w_i[jj]>0.01)) {
            res = beta_mean_inv - (1./b[jj]) ;
            if (fabs(res*w_i[jj])<resmax)          {;
            chi2 = chi2 + pow((res*w_i[jj]),2) ;
            icount= icount+1;
            sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
            sw =sw + w_i[jj]*w_i[jj] ;
                                               }
                                                                        }
                                      }
      quality = sqrt(sw) ;

      if (icount==0) chi2 = 1000.;
      if (icount>0) chi2 = chi2/(icount) ;
      if (icount>0) betachi=1./(sxw/sw);

   bxx = 100.;
   if ((chi2 < chi2cut)&&(quality>qualitycut)) bxx = betachi;
  //
  return(bxx);
};


////////////////////////////////////////////////////
////////////////////////////////////////////////////


/**
 * Fills a struct cToFLevel2 with values from a ToFLevel2 object (to put data into a F77 common).
 */
void ToFLevel2::GetLevel2Struct(cToFLevel2 *l2) const{

  for(Int_t i=0;i<6;i++)
    l2->tof_j_flag[i]=tof_j_flag[i];

  if(ToFTrk){ //ELENA
      l2->ntoftrk = ToFTrk->GetEntries();
      for(Int_t j=0;j<l2->ntoftrk;j++){
          l2->toftrkseqno[j]= ((ToFTrkVar*)ToFTrk->At(j))->trkseqno;
          l2->npmttdc[j]= ((ToFTrkVar*)ToFTrk->At(j))->npmttdc;
          for(Int_t i=0;i<l2->npmttdc[j];i++){
              l2->pmttdc[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->pmttdc.At(i);
              l2->tdcflag[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->tdcflag.At(i); // gf: 30 Nov 2006
          }
          for(Int_t i=0;i<13;i++)
              l2->beta[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->beta[i];
          
          l2->npmtadc[j]= ((ToFTrkVar*)ToFTrk->At(j))->npmtadc;
          for(Int_t i=0;i<l2->npmtadc[j];i++){
              l2->pmtadc[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->pmtadc.At(i);
              l2->adcflag[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->adcflag.At(i); // gf: 30 Nov 2006
              l2->dedx[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->dedx.At(i);
          }
          for(Int_t i=0;i<3;i++){
              l2->xtofpos[i][j]=((ToFTrkVar*)ToFTrk->At(j))->xtofpos[i];
              l2->ytofpos[i][j]=((ToFTrkVar*)ToFTrk->At(j))->ytofpos[i];
          }
          for(Int_t i=0;i<6;i++){
              l2->xtr_tof[i][j]=((ToFTrkVar*)ToFTrk->At(j))->xtr_tof[i];
              l2->ytr_tof[i][j]=((ToFTrkVar*)ToFTrk->At(j))->ytr_tof[i];
          }
      }
  } //ELENA
    
  if(PMT){ //ELENA
      l2->npmt = PMT->GetEntries();
      for(Int_t j=0;j<l2->npmt;j++){
          l2->pmt_id[j] = ((ToFPMT*)PMT->At(j))->pmt_id;
          l2->adc[j] =((ToFPMT*)PMT->At(j))->adc;
          l2->tdc_tw[j] =((ToFPMT*)PMT->At(j))->tdc_tw;
      }
  } //ELENA
}


//
// Reprocessing tool // Emiliano 08/04/07
// 
Int_t ToFLevel2::Process(TrkLevel2 *trk, TrigLevel2 *trg, GL_RUN *run, OrbitalInfo *orb, Bool_t force){
  // 
  // Copiare qui qualcosa di simile a calonuclei per evitare di riprocessare sempre tutto
  //


  //
  // structures to communicate with F77
  //
  extern struct ToFInput  tofinput_;
  extern struct ToFOutput tofoutput_;
  // 
  // DB connection
  //
  TString host;
  TString user;
  TString psw;
  const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
  const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
  const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
  if ( !pamdbhost ) pamdbhost = "";
  if ( !pamdbuser ) pamdbuser = "";
  if ( !pamdbpsw ) pamdbpsw = "";
  if ( strcmp(pamdbhost,"") ) host = pamdbhost;
  if ( strcmp(pamdbuser,"") ) user = pamdbuser;
  if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
  //
  //
  TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
  if ( !dbc->IsConnected() ) return 1;
  stringstream myquery;
  myquery.str("");
  myquery << "SET time_zone='+0:00'";
  dbc->Query(myquery.str().c_str());
  GL_PARAM *glparam = new GL_PARAM();
  glparam->Query_GL_PARAM(1,1,dbc); // parameters stored in DB in GL_PRAM table
  trk->LoadField(glparam->PATH+glparam->NAME);
  //
  Bool_t defcal = true;
  Int_t error=glparam->Query_GL_PARAM(run->RUNHEADER_TIME,201,dbc); // parameters stored in DB in GL_PRAM table
  if ( error<0 ) {
    return(1);
  };
  printf(" Reading ToF parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
  if ( (UInt_t)glparam->TO_TIME != (UInt_t)4294967295UL ) defcal = false;
  //
  Int_t nlen = (Int_t)(glparam->PATH+glparam->NAME).Length();
  rdtofcal((char *)(glparam->PATH+glparam->NAME).Data(),&nlen);
  //
  Int_t adc[4][12];
  Int_t tdc[4][12];
  Float_t tdcc[4][12];
  //
  // process tof data
  //
  for (Int_t hh=0; hh<12;hh++){
    for (Int_t kk=0; kk<4;kk++){
           adc[kk][hh] = 4095;
           tdc[kk][hh] = 4095;
           tdcc[kk][hh] = 4095.;
           tofinput_.adc[hh][kk] = 4095;
           tofinput_.tdc[hh][kk] = 4095;
    };
  };
  Int_t ntrkentry = 0;
  Int_t npmtentry = 0;
  Int_t gg = 0;
  Int_t hh = 0;
  Int_t adcf[48];
  memset(adcf, 0, 48*sizeof(Int_t));
  Int_t tdcf[48];
  memset(tdcf, 0, 48*sizeof(Int_t));
  for (Int_t pm=0; pm < this->ntrk() ; pm++){
     ToFTrkVar *ttf = this->GetToFTrkVar(pm);
     for ( Int_t nc=0; nc < ttf->npmttdc; nc++){
            if ( (ttf->tdcflag).At(nc) != 0 ) tdcf[(ttf->pmttdc).At(nc)] = 1;
     };
     for ( Int_t nc=0; nc < ttf->npmtadc; nc++){
            if ( (ttf->adcflag).At(nc) != 0 ) adcf[(ttf->pmtadc).At(nc)] = 1;
     };
  };
  //
  for (Int_t pm=0; pm < this->npmt() ; pm++){
     ToFPMT *pmt = this->GetToFPMT(pm);
     this->GetPMTIndex(pmt->pmt_id, gg, hh);
     if ( adcf[pmt->pmt_id] == 0 ){
             tofinput_.adc[gg][hh] = (int)pmt->adc;
             adc[hh][gg] = (int)pmt->adc;
     };
     if ( tdcf[pmt->pmt_id] == 0 ){
             tofinput_.tdc[gg][hh] = (int)pmt->tdc;
             tdc[hh][gg] = (int)pmt->tdc;
     };
     tdcc[hh][gg] = (float)pmt->tdc_tw;
     // Int_t pppid = this->GetPMTid(hh,gg);
     //      printf(" pm %i pmt_id %i pppid %i hh %i gg %i tdcc %f tdc %f adc %f \n",pm,pmt->pmt_id,pppid,hh,gg,pmt->tdc_tw,pmt->tdc,pmt->adc);
  };
  //
  Int_t unpackError = this->unpackError;
  //
  for (Int_t hh=0; hh<5;hh++){
     tofinput_.patterntrig[hh]=trg->patterntrig[hh];
  };
  //
  this->Clear();
  //
      Int_t pmt_id = 0;
      ToFPMT *t_pmt = new ToFPMT();
      if(!(this->PMT)) this->PMT = new TClonesArray("ToFPMT",12); //ELENA
      TClonesArray &tpmt = *this->PMT;
      ToFTrkVar *t_tof = new ToFTrkVar();
      if(!(this->ToFTrk)) this->ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
      TClonesArray &t = *this->ToFTrk;
      //
      //
      // Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract common and track-related  variables.
      //
      npmtentry = 0;
      //
      ntrkentry = 0;
      //
      // Calculate tracks informations from ToF alone
      //
      tofl2com();
      //
      memcpy(this->tof_j_flag,tofoutput_.tof_j_flag,6*sizeof(Int_t));
      //
      t_tof->trkseqno = -1;
      //
      // and now we must copy from the output structure to the level2 class:
      //
      t_tof->npmttdc = 0;
      //
      for (Int_t hh=0; hh<12;hh++){
        for (Int_t kk=0; kk<4;kk++){
          if ( tofoutput_.tofmask[hh][kk] != 0 ){
            pmt_id = this->GetPMTid(kk,hh);
            t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
            t_tof->tdcflag.AddAt(tofoutput_.tdcflagtof[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
            t_tof->npmttdc++;
          };
        };
      };
      for (Int_t kk=0; kk<13;kk++){
        t_tof->beta[kk] = tofoutput_.betatof_a[kk];
      }
      //
      t_tof->npmtadc = 0;
      for (Int_t hh=0; hh<12;hh++){
        for (Int_t kk=0; kk<4;kk++){
          if ( tofoutput_.adctof_c[hh][kk] < 1000 ){
            t_tof->dedx.AddAt(tofoutput_.adctof_c[hh][kk],t_tof->npmtadc);
            pmt_id = this->GetPMTid(kk,hh);
            t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
            t_tof->adcflag.AddAt(tofoutput_.adcflagtof[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
            t_tof->npmtadc++;
          };
        };
      };
      //
      memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
      memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
      memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
      memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
      //
      new(t[ntrkentry]) ToFTrkVar(*t_tof);
      ntrkentry++;
      t_tof->Clear();
      //
      //
      //
      t_pmt->Clear();
      //
      for (Int_t hh=0; hh<12;hh++){
        for (Int_t kk=0; kk<4;kk++){
         // new WM
          if ( tofoutput_.tdc_c[hh][kk] < 4095 || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){
//          if ( tdcc[kk][hh] < 4095. || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){
            //
            t_pmt->pmt_id = this->GetPMTid(kk,hh);
            t_pmt->tdc_tw = tofoutput_.tdc_c[hh][kk];
            t_pmt->adc = (Float_t)adc[kk][hh];
            t_pmt->tdc = (Float_t)tdc[kk][hh];
            //
            new(tpmt[npmtentry]) ToFPMT(*t_pmt);
            npmtentry++;
            t_pmt->Clear();
          };
        };
      };
      //
      // Calculate track-related variables
      //
      if ( trk->ntrk() > 0 ){
        //
        // We have at least one track
        //
        //
        // Run over tracks
        //
        for(Int_t nt=0; nt < trk->ntrk(); nt++){
          //
          TrkTrack *ptt = trk->GetStoredTrack(nt);
          //
          // Copy the alpha vector in the input structure
          //
          for (Int_t e = 0; e < 5 ; e++){
            tofinput_.al_pp[e] = ptt->al[e];
          };
          //
          // Get tracker related variables for this track
          //
          toftrk();
          //
          // Copy values in the class from the structure (we need to use a temporary class to store variables).
          //
          t_tof->npmttdc = 0;
          for (Int_t hh=0; hh<12;hh++){
            for (Int_t kk=0; kk<4;kk++){
              if ( tofoutput_.tofmask[hh][kk] != 0 ){
                pmt_id = this->GetPMTid(kk,hh);
                t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
                t_tof->tdcflag.AddAt(tofoutput_.tdcflag[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
                t_tof->npmttdc++;
              };
            };
          };
          for (Int_t kk=0; kk<13;kk++){
            t_tof->beta[kk] = tofoutput_.beta_a[kk];
          };
          //
          t_tof->npmtadc = 0;
          for (Int_t hh=0; hh<12;hh++){
            for (Int_t kk=0; kk<4;kk++){
              if ( tofoutput_.adc_c[hh][kk] < 1000 ){
                t_tof->dedx.AddAt(tofoutput_.adc_c[hh][kk],t_tof->npmtadc);
                pmt_id = this->GetPMTid(kk,hh);
                t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
                t_tof->adcflag.AddAt(tofoutput_.adcflag[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
                t_tof->npmtadc++;
              };
            };
          };
          //
          memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
          memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
          memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
          memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
          //
          // Store the tracker track number in order to be sure to have shyncronized data during analysis
          //
          t_tof->trkseqno = nt;
          //
          // create a new object for this event with track-related variables
          //
          new(t[ntrkentry]) ToFTrkVar(*t_tof);
          ntrkentry++;
          t_tof->Clear();
          //
        }; // loop on all the tracks
      //
      this->unpackError = unpackError;
      if ( defcal ){
        this->default_calib = 1;
      } else {
        this->default_calib = 0;
      };
 };


  return(0);
}