TrackerLevel2/src/ExtTrkingAlg.cpp

/**
 * \file ExtTrkingAlg.cpp
 * \author Elena Vannuccini
 **/
#include <ExtTrkingAlg.h>
#include <iostream>
#include <math.h>

using namespace std;
/**
 * Constructor
 *
 * @param id Algorythm identifier
 *
 * id = 0 The algorythm looks for tracks with >= 3x2y hits on the tracker plane. 
 *        It fills a TClonesArray of TrkTrack objects. 
 * id = 1 The first plane (y) of the calorimeter is also used. 
 *        It fills a TClonesArray of ExtTrack objects. 
 * id = 2 The first two planes (y+x) of the calorimeter are also used. 
 *        It fills a TClonesArray of ExtTrack objects. 
 *
 */
ExtTrkingAlg::ExtTrkingAlg(Int_t id){

  _whichAlg = id; 

  SetDebug(false);

   
  if(id == 0){
    if(_debug){
      cout << "ExtTrkingAlg::ExtTrkingAlg("<<id<<")"<<endl;
      cout << "Creating array of TrkTrack objects "<<endl;
      cout << "WARNING!!! tracking not accurate!! bug not fixed yet... "<<endl;
    }
    _trkArray = new TClonesArray("TrkTrack");
    //
    _sel_nClstrMAX  = 20. ; // selection parameter: maximum number of cluster
    _sel_nPlaneXMIN = 3.; // selection parameter: minimum number of hit x-views
    _sel_nPlaneYMIN = 2.; // selection parameter: minimum number of hit y-views
    //
    _alg_nClFixX   = 3.;     // algorythm parameter: n.hits required on X view
    _alg_nClFixY   = 2.;     // algorythm parameter:n.hits required on Y view
    _alg_nTrackMAX = 400.;     // algorythm parameter: maximum num. of track candidates
    _alg_nViewCal= 0; // ignored
    //    _alg_nViewCalFix= 0; // ignored

    NEXTVIEWS  = 12 +_alg_nViewCal;
    NEXTPLANES = 6  +_alg_nViewCal;
    _extTrack = NULL;//ignored
    _zMech = NULL;//ignored
    

  // }else if (id == 1 ){
  }else if ( (id >= 100 && id < 144) || 
             (id >= 200 && id < 244) || 
             false){ //add n calo planes
    if(_debug){
      cout << "ExtTrkingAlg::ExtTrkingAlg("<<id<<")"<<endl;
      cout << "Creating array of ExtTrack objects "<<endl;
    }
    _trkArray = new TClonesArray("ExtTrack");
    //
    _alg_nClFixX   = 3.;     // algorythm parameter: n.hits required on X view
    _alg_nClFixY   = 2.;     // algorythm parameter:n.hits required on Y view
    _alg_nTrackMAX = 1000.;     // algorythm parameter: maximum num. of track candidates
    // _alg_nViewCal = 1;        // algorythm parameter: n.calorimeter planes included
    //_alg_nViewCal = 1;        // algorythm parameter: n.calorimeter planes included
    
    if((int)(id/100)==1)_alg_nViewCal = id - 100;  // algorythm parameter: n.calorimeter planes included
    if((int)(id/200)==1)_alg_nViewCal = id - 200;  // algorythm parameter: n.calorimeter planes included
 
    //
    _sel_nClstrMAX  = 15.;//10. ; // selection parameter: maximum number of cluster
    _sel_nPlaneXMIN = (_alg_nViewCal ? 2 : 3);   // selection parameter: minimum number of hit x-views
    _sel_nPlaneYMIN = 2 - (int)(_alg_nViewCal+1)/2 ;   // selection parameter: minimum number of hit y-views

    NEXTVIEWS  = 12 +_alg_nViewCal;
    NEXTPLANES = 6  +_alg_nViewCal;
    _extTrack = new ExtTrack(NEXTPLANES);
    _zMech = new float[NEXTPLANES];
    
    // -----------------------------------
    // set Z-coord of tracker planes
    // -----------------------------------  
    _zMech[0]=ZTRK1;
    _zMech[1]=ZTRK2;
    _zMech[2]=ZTRK3;
    _zMech[3]=ZTRK4;
    _zMech[4]=ZTRK5;
    _zMech[5]=ZTRK6;
    // _extTrack->SetZ(0,ZTRK1);
    // _extTrack->SetZ(1,ZTRK2);
    // _extTrack->SetZ(2,ZTRK3);
    // _extTrack->SetZ(3,ZTRK4);
    // _extTrack->SetZ(4,ZTRK5);
    // _extTrack->SetZ(5,ZTRK6);
    
    _caloStripRoto.clear();
    for(int iv=0; iv<_alg_nViewCal; iv++){ //loop over calorimater tracking planes 
      for(int is=0; is<9; is++){ //loop over calorimeter sensors
          CaloStripRoto s(1-iv%2, (int)(iv/2), is, false);//flight 0-order align         
          s.ResetAligParams();//reset 1-order align parameters   
          _caloStripRoto.push_back(s) ;
      }
    }
    //------------------------------------------
    // read alignment parameters
    //------------------------------------------
    if(_debug) cout << "Reading calorimeter alignment parameters"<<endl;
    const char *pamca = gSystem->Getenv("PAM_CALIB");
    TString filep = "/trk-param/align_param_calo-0/CaloAlignParams.txt";
    filep.Prepend(pamca);
    std::fstream fs;
    fs.open (filep.Data(), std::fstream::in );
    do{
        int view,plane,sensor;
        float par[5],err[5];
        fs>>view >> plane>> sensor;
        for(int i=0; i<5; i++)fs>>par[i]>>err[i];
        
        
        if(plane*2+1-view<_alg_nViewCal){


            int index = (plane*2+1-view)*9 + sensor;
            if(_debug){
              cout<<index<<" >> ";
              cout<<view << " "<<plane<<" "<<sensor<<"          ";
              for(int i=0; i<5; i++)cout<<setw(10)<<par[i]<<" "<<setw(10)<<err[i]<<" ";
              cout<<endl;
            }

            _caloStripRoto[index].SetAligParams__fPitch(par[0]);
            _caloStripRoto[index].SetAligParams__shift(par+1);
            _caloStripRoto[index].SetAligParams__alpha(par[4]);
//          _caloStripRoto[index].SetSiSensor(sensor);
//          _caloStripRoto[index].SetView(view);
//          _caloStripRoto[index].SetPlane(plane);

        }
    }while(!fs.eof());
    fs.close();

    // -----------------------------------
    // set Z-coord of calorimeter planes
    // -----------------------------------      
    for(int plane=0; plane<22; plane++){
      for(int view=1; view>=0; view--){
        int strip = 0;
        CaloStrip st = CaloStrip(false);//flight align
        st.Set(view,plane,strip);       
        int index = 2*plane + !view;
        if(index >= _alg_nViewCal)break;
        // _extTrack->SetZ(6+index,st.GetZ());
        _zMech[6+index]=st.GetZ();
      }
    }

    for(int ip =0 ; ip<_extTrack->nplanes; ip++ )_extTrack->SetZ(ip,_zMech[ip]);

    if ( _debug ){
      cout <<" Extended-track Z-coordinates: "<<endl;
      for(int ip =0 ; ip<_extTrack->nplanes; ip++ )cout <<  _extTrack->zm[ip]<<endl;
    }
  }else{
    if ( _debug ){
      cout << "ExtTrkingAlg(Int_t id) -- algorythm id="<<id<<" not valid "<<endl;
      cout << "--> Track array not created "<<endl;
    }
    _sel_nClstrMAX  = 0.; 
    _sel_nPlaneXMIN = 0.;   
    _sel_nPlaneYMIN = 0.;  
    //
    _alg_nClFixX   = 0.;   
    _alg_nClFixY   = 0.;  
    _alg_nTrackMAX = 0.;   
    _alg_nViewCal = 0;     
    //    _alg_nViewCalFix= 0; 

    NEXTVIEWS  = 12 +_alg_nViewCal;
    NEXTPLANES = 6  +_alg_nViewCal;
    _extTrack = NULL;//ignored

  }

  CaloStrip st = CaloStrip(false);
  for(int view=0; view<2; view++){
    for(int plane=0; plane<22; plane++){
      for(int strip=0; strip<96; strip++){
        st.Set(view,plane,strip);
        _caloCoord[view][plane][strip]=(view ? st.GetY() : st.GetX());
        _caloZ[plane*2+(view?0:1)]=st.GetZ();
      }
    }
  }
  _caloTj = new Trajectory(44,_caloZ);
//  _caloTj->Dump();

  int    ngf = TrkParams::nGF;
  float *zgf = TrkParams::zGF;
  _tgf = new Trajectory(ngf,zgf);
//  _tgf->Dump();
 
  Clear();
  
};

/**
 * Reset the algorythm flags and output
 */
// void ExtTrkingAlg::Reset(){

//   if(_trkArray)_trkArray->Clear("C");

// };

/**
 * Clear the event and reset the algorythm
 */
void ExtTrkingAlg::Clear(Option_t* option){

  SetTrkLevel1();
  SetTrkLevel2();

  SetToFLevel2();

  SetCaloLevel1();
  SetCaloLevel2();

  if(_trkArray)_trkArray->Clear("C");

};  
/**
 * DElete all
 */

void ExtTrkingAlg::Delete(){ 
  if(_trkArray) delete _trkArray; 
  if(_caloTj) delete _caloTj;
  if(_tgf) delete _tgf;
  if(_zMech)delete [] _zMech;
}


/**
 * Set selection parameters
 */
void ExtTrkingAlg::SetSelectionParams(double* par){

  _sel_nClstrMAX  = (Int_t)par[0]; // selection parameter: maximum number of cluster
  _sel_nPlaneXMIN = (Int_t)par[1]; // selection parameter: minimum number of hit x-views
  _sel_nPlaneYMIN = (Int_t)par[2]; // selection parameter: minimum number of hit y-views

};
/**
 * Set algorythm parameters
 */
void ExtTrkingAlg::SetAlgorythmParams(double* par){

  _alg_nClFixX   = (Int_t)par[0];     // algorythm parameter: n.hits required on X view
  _alg_nClFixY   = (Int_t)par[1];     // algorythm parameter:n.hits required on X view
  _alg_nTrackMAX = (Int_t)par[2];     // algorythm parameter: maximum num. of track candidates
  _alg_nViewCal=  (Int_t)par[3];     // algorythm parameter: n. calorimeter planes included 
  //  _alg_nViewCalFix=  (Int_t)par[4];     // 
 
};

/**
 * Get the track array
 */
TClonesArray* ExtTrkingAlg::GetTrackArray(Bool_t reset){
  
  if( reset && _trkArray) _trkArray->Clear("C");
  //  if( !_procDone )ProcessEvent();

  return _trkArray;

};
/**
 * Pre-select the event
 */
Bool_t ExtTrkingAlg::CheckEvent(){

  if(!_trk_l2)return true;

  TrkLevel2 *trkl2 = _trk_l2;//event->GetTrkLevel2();
  int nClstrMAX  = _sel_nClstrMAX; //maximum number of cluster
  int nPlaneXMIN = _sel_nPlaneXMIN; 
  int nPlaneYMIN = _sel_nPlaneYMIN;
  /////////////////////////////////////////////////////////////////////
  /// count number of hit planes
  /////////////////////////////////////////////////////////////////////
  int nHitX[] = {0,0,0,0,0,0};
  int nHitY[] = {0,0,0,0,0,0};
  for(int ix=0; ix<trkl2->nclsx(); ix++)nHitX[trkl2->GetSingletX(ix)->plane-1]++;    
  for(int iy=0; iy<trkl2->nclsy(); iy++)nHitY[trkl2->GetSingletY(iy)->plane-1]++;
  /////////////////////////////////////////////////////////////////////
  /// set minimum condition 3x+2y
  /////////////////////////////////////////////////////////////////////
  int nPlaneX=0;
  for(int ip=0; ip<6; ip++)if(nHitX[ip])nPlaneX++;      
  int nPlaneY=0;
  for(int ip=0; ip<6; ip++)if(nHitY[ip])nPlaneY++;      
  /////////////////////////////////////////////////////////////////////
  /// dump selection
  /////////////////////////////////////////////////////////////////////
  if(_debug){
    cout << " n.tracks "<<trkl2->GetNTracks()<<endl;
    cout << " n.x-singles "<<trkl2->nclsx()<<endl;
    cout << " n.y-singles "<<trkl2->nclsy()<<endl;
    cout << " n.x-planes with singles "<<nPlaneX<<endl;
    cout << " n.y-planes with singles "<<nPlaneY<<endl;
  }
  ////////////////////////////////////////////////////////////////////////
  // 
  // check condition for retracking
  //
  ////////////////////////////////////////////////////////////////////////
  if(trkl2->GetNTracks()==0 && nPlaneX<nPlaneXMIN) return false; 
  if(trkl2->GetNTracks()==0 && nPlaneY<nPlaneYMIN) return false; 
  if(trkl2->nclsy()+trkl2->nclsx() > nClstrMAX) return false; 

  return true;

};


/**
 * Process the event
 */
void ExtTrkingAlg::ProcessEvent0(Bool_t force){

  if(_debug)cout << " |---------------------------------------------------| "<<endl;
  if(_debug)cout << " void ExtTrkingAlg::ProcessEvent0("<<force<<") "<<endl;

  
  if(_debug && !_trk_l1)cout << " Missing TrkLevel1 object --- ciao! "<<endl;
  if(!_trk_l1)return;


  //===========================
  // pre-selection  
  //===========================
  
  Bool_t RETRACK = CheckEvent();

  if(_debug)cout <<" Pre-selection result: "<<RETRACK<<endl;
  if( !force && !RETRACK )return;
  

  //===========================
  // tracking
  //===========================

  // TClonesArray &tracks = *_trkArray;
  // tracks.Clear();


  multimap<int,TrkTrack> trackCandidates[12]; // map: last filled view vs TrkTrack
  int mapIndex = 0;
  multimap<int,TrkTrack>::iterator trkBegin;
  multimap<int,TrkTrack>::iterator trkEnd;
  int iCand;

  // ------------------------------------------------
  // fill a map VIEW(0-11)-CLUSTER
  // ------------------------------------------------
  multimap<int,int> clusterMap;
  if(_debug)cout<<"-------- CLUSTERs --------"<<endl;
    for(int icl=0; icl<_trk_l1->nclstr(); icl++){
    clusterMap.insert( pair<int,int>( _trk_l1->GetCluster(icl)->view-1, icl ));
    if(_debug)cout << icl << " - view "<<_trk_l1->GetCluster(icl)->view<<" ladder "<<_trk_l1->GetCluster(icl)->GetLadder()<<" strip "<<_trk_l1->GetCluster(icl)->maxs<<endl;
  }
  // ------------------------------------------------
  // define iterators over clusters on the same view
  // ------------------------------------------------
  pair <multimap<int,int>::iterator, multimap<int,int>::iterator> ret[12];
  multimap<int,int>::iterator it[12];
  for(int iv=0; iv<12; iv++)ret[iv] = clusterMap.equal_range(iv);
  
  //====================================================================
  // PARAMETERS
  //====================================================================
  // ------------------------------------------------
  // number of points per candidate
  // ------------------------------------------------
  int nClFixX =  _alg_nClFixX; // n.hits required on X view
  int nClFixY = _alg_nClFixY; // n.hits required on X view
  int nTrackMAX =  _alg_nTrackMAX; // 
  //==================================================
  // ------------------------------------------------
  // start: one candidate for each cluster
  // ------------------------------------------------
  //==================================================
  if(_debug)cout<<"-------- MINIMAL --------"<<endl;
  for(int iv=0; iv<12; iv++){                                      //loop over all views
    
    if( !(iv%2) && (int)(0.5*iv) >  6-nClFixX)continue;
    if(  (iv%2) && (int)(0.5*iv) >  6-nClFixY)continue;
    
    for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view
      for (int is=0; is<2; is++){                                  //loop over sensors
        int icl = it[iv]->second;
        TrkCluster *cluster = _trk_l1->GetCluster(icl);
        int ladder = cluster->GetLadder()-1;
        int plane = cluster->GetPlane()-1;
        bool isY = (cluster->view)%2;
        bool isX = !isY;
        //
        // if(_debug && isX )cout << "  ^^^^^^^ X plane "<<plane<<" ladder "<<ladder<<" strip "<<cluster->maxs <<endl;
        // if(_debug && isY )cout << "  ^^^^^^^ Y plane "<<plane<<" ladder "<<ladder<<" strip "<<cluster->maxs <<endl;
        
        TrkTrack t_track = TrkTrack();
        
        if( isX )t_track.SetXGood(plane,icl+1,ladder,is);       
        if( isY )t_track.SetYGood(plane,icl+1,ladder,is);
        
        trackCandidates[mapIndex].insert(pair<int,TrkTrack>(iv, TrkTrack(t_track)));
        
      }
    }
  };
  
  //==============================================================
  // -------------------------------------------------------------
  // next: add views to each candidate, up to nClFixX+nClFixY
  // -------------------------------------------------------------
  //==============================================================
  for( ;mapIndex < nClFixX+nClFixY-1 ;++mapIndex){

    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    if(_debug)cout<<"MAP "<<mapIndex<<" -----------"<<endl;

    for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){

      TrkTrack trackCand = (*trkIter).second;
      int lastView = (*trkIter).first; //last inserted view
      int lastPlane = (int)(0.5*lastView);
      int lastLadder = trackCand.GetLadder(lastPlane);
      int lastSensor = trackCand.GetSensor(lastPlane);

      // if(_debug)cout<<"MAP "<<mapIndex<<" last v "<<lastView<<" NX "<<trackCand.GetNX()<<" NY "<<trackCand.GetNY()<<endl;

      for(int iv=lastView+1; iv<12; iv++){                                      //loop over next views
                     
        if(  (iv%2) && (int)(0.5*iv) >  6-nClFixX+trackCand.GetNX())continue; //not enough x-view left
        if( !(iv%2) && (int)(0.5*iv) >  6-nClFixY+trackCand.GetNY())continue; //not enough y-view left

        if(  (iv%2) && trackCand.GetNX()==nClFixX )continue;//ok?
        if( !(iv%2) && trackCand.GetNY()==nClFixY )continue;//ok?

        for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view

          int icl = it[iv]->second;
          TrkCluster *cluster = _trk_l1->GetCluster(icl);
          int ladder = cluster->GetLadder()-1;//(_trk_l1->GetCluster(icl)->maxs-1)%1024 + 1;
          int plane = cluster->GetPlane()-1;
          bool isY = (cluster->view)%2;
          bool isX = !isY;
                
          if( plane==lastPlane && ladder!=lastLadder)continue;
                

          for (int is=0; is<2; is++){                                  //loop over sensors
                  
            if( plane==lastPlane && is!=lastSensor)continue;

            // if(_debug)cout<<"MAP "<<mapIndex<<" last v "<<lastView;
            // if(_debug)cout<<" +cl "<<icl<<" s"<<is<<endl;

            TrkTrack t_track = TrkTrack();
            trackCand.Copy(t_track);                             //copy the candidate parameters...
                  
            if( isX )t_track.SetXGood(plane,icl+1,ladder,is);         //...add a point...
            if( isY )t_track.SetYGood(plane,icl+1,ladder,is);
                  
            trackCandidates[mapIndex+1].insert(pair<int,TrkTrack>(iv, TrkTrack(t_track))); //...and store a new candidate
          };//end loop over sensors
        };//end loop over clusters on the view iv
      };//endl loop over views
    }//end loop over candidates
  }//end iterations
  

  if( trackCandidates[mapIndex].size() > (uint)nTrackMAX ){ // EM,  warning: comparison between signed and unsigned integer expressions [-Wsign-compare]
    if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<" > "<<nTrackMAX<<endl;
    return;//to many candidates
  }
  if(_debug){
    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
      TrkTrack &trackCand = (*trkIter).second;
      cout << " >> ";
      cout << " X: ";
      for(int ip=0; ip<6; ip++)
        if(trackCand.GetClusterX_ID(ip)>=0)cout << trackCand.GetClusterX_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";
      cout << " Y: ";
      for(int ip=0; ip<6; ip++)
        if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";     
      cout << endl;
    }
  }
  if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;

  //==============================================================
  // -------------------------------------------------------------
  // first selection of track candidates, based on chi2
  // -------------------------------------------------------------
  //==============================================================
  if(_debug)cout<<"-------- CHI2 --------"<<endl;
  trkBegin = trackCandidates[mapIndex].begin();
  trkEnd   = trackCandidates[mapIndex].end();
  mapIndex++;
  iCand = 0;
  for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
    int lastView = (*trkIter).first;
    TrkTrack &trackCand = (*trkIter).second;      
    int ifail=0;
    trackCand.FitReset();
    trackCand.Fit(0.,ifail,0,1);//re-evaluate cluster position from level1
    if(ifail!=0)trackCand.FitReset();

    // -----------------------------------
    // first selection of track candidates
    // -----------------------------------
    if(TMath::IsNaN(trackCand.chi2))continue;
    if(trackCand.chi2 <= 0. || trackCand.chi2 > 1.e7)continue; 
    //    if(trackCand.chi2 < 0. || trackCand.chi2 > 10. )continue; 
    iCand++;    
          
    if(_debug){     
     cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
      cout <<" | al: ";
      for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
      cout <<" | ";
      for(int ip=0; ip<6; ip++)
        if(trackCand.GetClusterX_ID(ip)>=0)cout << trackCand.GetClusterX_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";
      cout << " Y: ";
      for(int ip=0; ip<6; ip++)
        if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";       cout << endl;
    }

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // evaluated coordinates (to define GF)
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    int    ngf = TrkParams::nGF;
    float *zgf = TrkParams::zGF;
    Trajectory tgf = Trajectory(ngf,zgf);
    tgf.DoTrack(trackCand.al);//<<<< integrate the trajectory
    for(int ip=0; ip<ngf; ip++){
      trackCand.xGF[ip] = tgf.x[ip];
      trackCand.yGF[ip] = tgf.y[ip];
    }
        
    trackCandidates[mapIndex].insert(pair<int,TrkTrack>(lastView, TrkTrack(trackCand)));

  };
  if(trackCandidates[mapIndex].size()==0)return;//no good candidates

  if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
  //=================================================================
  // ----------------------------------------------------------------
  // CALORIMETER
  // ----------------------------------------------------------------
  //=================================================================
  if(_cal_l1){
    if(_debug)cout<<"-------- CALORIMETER --------"<<endl;
    // -----------------------------------
    // evaluate calorimeter variables
    // -----------------------------------
    // float caloCoord[2][22][96];
    // float caloZ[44];
    // for(int view=0; view<2; view++){
    //   for(int plane=0; plane<22; plane++){
    //  for(int strip=0; strip<96; strip++){
    //    CaloStrip st = CaloStrip(false);
    //    st.Set(view,plane,strip);
    //    caloCoord[view][plane][strip]=(view ? st.GetY() : st.GetX());
    //    caloZ[plane*2+(view?0:1)]=st.GetZ();
    //  }
    //   }
    // }
    // Trajectory caloTj = Trajectory(44,caloZ);
    //  for(int iz=0; iz<44;iz++)cout<<caloZ[iz]<<endl;
    //test calo
    //    float caloEvent[2][22][96];
    CaloLevel1* l1=_cal_l1;//event->GetCaloLevel1();
    //    if(!l1)return;
    vector<float> caloChi2;
    vector<int> caloHit; int caloHitMAX=0;
    vector<float> caloMip;
    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){        
      TrkTrack trackCand = ((*trkIter).second);   
            
      trackCand.DoTrack(_caloTj);//integrate the trajectory through the calorimeter
            
      float chi2 = 0.;
      int  nhit = 0;
      float cmip = 0.;
      for(Int_t ih=0; ih<l1->istrip; ih++){
        Int_t view = -1;
        Int_t plane = -1;
        Int_t strip = -1;
        Float_t mip = l1->DecodeEstrip(ih,view,plane,strip);
        if(view<0)continue;
        if(plane<0)continue;
        if(strip<0)continue;
        if(mip<=0)continue;
        float dxy = _caloCoord[view][plane][strip]-(view?_caloTj->y[plane*2+(view?0:1)]:_caloTj->x[plane*2+(view?0:1)]);
        if(plane<11) chi2+= mip * TMath::Power(dxy,2.);
        //            caloZ[plane*2+(view?0:1)]=st.GetZ();
        if( fabs(dxy) < 2. && plane<11)nhit++;
        if( fabs(dxy) < 2. && plane<11)cmip+=mip;
      };
      if(l1->istrip>0)chi2 = chi2/l1->istrip;
            
      caloHit.push_back(nhit); 
      if(nhit>caloHitMAX)caloHitMAX=nhit;

      if(_debug){           
        cout << " CALO ";
        //      cout << " chi2 "<< chi2 <<" ";
        cout << " nhit (< 11th plane) "<< nhit <<" ";
        //      cout << " mip  "<< cmip <<" ";
        cout <<endl;
      };
    }
          
    //=================================================================
    // ----------------------------------------------------------------
    // selection of track candidates, based on n.calorimeter hits
    // ----------------------------------------------------------------
    //=================================================================
    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    mapIndex++;
    iCand = 0;
    for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
      int lastView = (*trkIter).first;
      TrkTrack &trackCand = (*trkIter).second;
      if(caloHit[iCand]==caloHitMAX)trackCandidates[mapIndex].insert(pair<int,TrkTrack>(lastView, TrkTrack(trackCand)));
      iCand++;
    }
    if(trackCandidates[mapIndex].size()==0)return;//no good candidates
    if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;

  }

  //=================================================================
  // ----------------------------------------------------------------
  // TOF
  // ----------------------------------------------------------------
  //=================================================================
  if(_tof_l2){
    if(_debug)cout<<"-------- TOF --------"<<endl;
    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    ToFLevel2* tl2=_tof_l2;//event->GetToFLevel2(); 
    //   if(!tl2)return;
    iCand =0;
    vector<int> tofHit; int tofHitMAX=0; 
    for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
      //            int lastView = (*trkIter).first;
      TrkTrack &trackCand = (*trkIter).second;
            
      int nhit =0;
      for(int iToF=0; iToF<6; iToF++){
              
        int iGF = iToF;
        if(iToF>3)iGF = iToF+8;
        int iPaddle = tl2->GetPaddleIdOfTrack( trackCand.xGF[iGF], trackCand.yGF[iGF], iToF);
        if( tl2->HitPaddle(iToF,iPaddle) )nhit++;
              
      }
      iCand++;
            
      tofHit.push_back(nhit); 
      if(nhit>tofHitMAX)tofHitMAX=nhit;

      if(_debug){
        cout << " TOF nhit "<< nhit <<" ";
        cout <<endl;
      }
            
    }
    //=================================================================
    // ----------------------------------------------------------------
    // selection of track candidates, based on n.tof hits
    // ----------------------------------------------------------------
    //=================================================================
    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    mapIndex++;
    iCand = 0;
    for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
      int lastView = (*trkIter).first;
      TrkTrack &trackCand = (*trkIter).second;
      if(tofHit[iCand]==tofHitMAX)trackCandidates[mapIndex].insert(pair<int,TrkTrack>(lastView, TrkTrack(trackCand)));
      iCand++;
    }
    if(trackCandidates[mapIndex].size()==0)return;//no good candidates
    if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;

  }


  //=================================================================
  // ----------------------------------------------------------------
  // TRACK MERGING
  // ----------------------------------------------------------------
  //=================================================================
  if(_debug)cout<<"-------- MERGING --------"<<endl;
  trkBegin = trackCandidates[mapIndex].begin();
  trkEnd   = trackCandidates[mapIndex].end();
  mapIndex++;
  float chi2Min=1e10;
  iCand = 0;
  for( multimap<int,TrkTrack>::iterator trkIter1 = trkBegin; trkIter1!=trkEnd; ++trkIter1){ //loop over tracks (1)
    
    TrkTrack &trackCand1 = (*trkIter1).second;// get candidate 1...
            
    TrkTrack t_track = TrkTrack();// ...create a new TrkTrack instance...
    trackCand1.Copy(t_track);     // ...and make a copy
                   
    bool HBM = false; //has been merged

    for( multimap<int,TrkTrack>::iterator trkIter2 = trkBegin; trkIter2!=trkEnd; ++trkIter2){//loop over tracks (2)
              
      TrkTrack &trackCand2 = (*trkIter2).second; // get candidate 2...

      // if(_debug){
      //        cout <<" X: ";
      //        for(int ip=0; ip<6; ip++)
      //          if(t_track.GetClusterX_ID(ip)>=0)cout << t_track.GetClusterX_ID(ip)<<":"<<t_track.GetSensor(ip)<<"|";
      //        cout << " Y: ";
      //        for(int ip=0; ip<6; ip++)
      //          if(t_track.GetClusterY_ID(ip)>=0)cout << t_track.GetClusterY_ID(ip)<<":"<<t_track.GetSensor(ip)<<"|";         
      //        cout << " + ";
      //        cout <<" X: ";
      //        for(int ip=0; ip<6; ip++)
      //          if(trackCand2.GetClusterX_ID(ip)>=0)cout << trackCand2.GetClusterX_ID(ip)<<":"<<trackCand2.GetSensor(ip)<<"|";
      //        cout << " Y: ";
      //        for(int ip=0; ip<6; ip++)
      //          if(trackCand2.GetClusterY_ID(ip)>=0)cout << trackCand2.GetClusterY_ID(ip)<<":"<<trackCand2.GetSensor(ip)<<"|";        
      //        cout << " -----> ";
      // }

      bool CBM = true;
      for(int ip=0; ip<6; ip++){ //loop over planes
        if( t_track.XGood(ip)  &&  
            trackCand2.XGood(ip) &&  
            ( t_track.GetClusterX_ID(ip)!=trackCand2.GetClusterX_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
            true )CBM=false; // if both have a x-hit and it is not the same ---> they cannot be merged
        if(!CBM)break;       
        if( t_track.YGood(ip)  && 
            trackCand2.YGood(ip) && 
            ( t_track.GetClusterY_ID(ip)!=trackCand2.GetClusterY_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
            true )CBM=false; // if both have a y-hit and it is not the same ---> they cannot be merged
        if(!CBM)break;        
      }
      //      if(_debug)cout << " can be merged ? "<<CBM<<endl;
      if(!CBM)continue; //go to the next candidate
      ///////////
      //merging//
      ///////////      
      for(int ip=0; ip<6; ip++){              
        if( !(t_track.XGood(ip)) &&  trackCand2.XGood(ip) )
          t_track.SetXGood(ip,trackCand2.GetClusterX_ID(ip)+1,trackCand2.GetLadder(ip),trackCand2.GetSensor(ip));       
        if( !(t_track.YGood(ip)) &&  trackCand2.YGood(ip) ) 
          t_track.SetYGood(ip,trackCand2.GetClusterY_ID(ip)+1,trackCand2.GetLadder(ip),trackCand2.GetSensor(ip));
      };         
      HBM = true;
    }//end loop over tracks (2)


    // cout << "<<  ";
    // for(int ip=0; ip<6; ip++) cout <<setw(4)<< t_track.GetClusterX_ID(ip);
    // for(int ip=0; ip<6; ip++) cout <<setw(4)<< t_track.GetClusterY_ID(ip);
    // cout << endl;
    // cout<< " TRK chi2 "<<setw(10)<<t_track.chi2;
    // cout <<" | al ";
    // for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
    // cout << endl;
 
    ///////////////////////////////////////////////////////// FIT START
    // int ifail=0;
    // //           TrkParams::SetDebugMode();
    // t_track.FitReset();       //angles set to zero
    // t_track.Fit(0.,ifail,0,1);//re-evaluate cluster position from level1
    // if( TMath::IsNaN(t_track.chi2) )continue;
    // if( ifail != 0 )continue;
    // t_track.Fit(0.,ifail,0,1);//evaluate again taking into account the track angle
    // if( TMath::IsNaN(t_track.chi2) )continue;
    // if( ifail != 0 )continue;
    ///////////////////////////////////////////////////////// FIT FINISH
    // cout << "Merged: "<<endl;
    // cout << "X "; for(int ip=0; ip<6; ip++)cout << t_track.XGood(ip); cout << endl;
    // cout << "Y "; for(int ip=0; ip<6; ip++)cout << t_track.YGood(ip); cout << endl;
   //
    //check if track has been already inserted
    //
    multimap<int,TrkTrack>::iterator trkB = trackCandidates[mapIndex].begin();
    multimap<int,TrkTrack>::iterator trkE = trackCandidates[mapIndex].end();
    //
    bool ADD = true;
    for( multimap<int,TrkTrack>::iterator trkI = trkB; trkI!=trkE; ++trkI){
      TrkTrack &trackC = (*trkI).second;
      bool SAME=true;
      for(int ip=0; ip<6; ip++)
        if( t_track.GetClusterX_ID(ip) != trackC.GetClusterX_ID(ip) || 
            t_track.GetClusterY_ID(ip) != trackC.GetClusterY_ID(ip) || 
            t_track.GetSensor(ip)      != trackC.GetSensor(ip)      ||
            false)SAME=false;
      if(SAME){
        ADD=false;
        break;
      }
    }
    //    if(_debug)cout <<" ADD ? "<<ADD<<endl;
    if(ADD){

      //      cout << "-----> INSERT  "<<endl;
      if(_debug){
        cout << " >> HBM "<<HBM;
        cout << " X:";for(int ip=0; ip<6; ip++)cout << t_track.XGood(ip);
        cout << " Y:";for(int ip=0; ip<6; ip++)cout << t_track.YGood(ip);       
        cout << " S:";for(int ip=0; ip<6; ip++)cout << t_track.GetSensor(ip);   
      }

      ///////////////////////////////////////////////////////// FIT START
      int ifail=0;      
      t_track.FitReset();       //angles set to zero
      t_track.Fit(0.,ifail,0,1);//re-evaluate cluster position from level1
      if( TMath::IsNaN(t_track.chi2) )continue;
      if( ifail != 0 )continue;
      t_track.Fit(0.,ifail,0,1);//evaluate again taking into account the track angle
      if( TMath::IsNaN(t_track.chi2) )continue;
      if( ifail != 0 )continue;
      ///////////////////////////////////////////////////////// FIT FINISH

      if(_debug){
        cout << " >> al: ";for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
        cout << " chi2: "<<setw(10)<< t_track.chi2;
        cout << endl;
      }
      
      // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      // evaluated coordinates (to define GF)
      // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      int    ngf = TrkParams::nGF;
      float *zgf = TrkParams::zGF;
      Trajectory tgf = Trajectory(ngf,zgf);
      tgf.DoTrack(t_track.al);//<<<< integrate the trajectory
      for(int ip=0; ip<ngf; ip++){
        t_track.xGF[ip] = tgf.x[ip];
        t_track.yGF[ip] = tgf.y[ip];
      }
      // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      // evaluated other track info
      // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~      
      for(int ip=0; ip<6; ip++){
        float factor = 0.;
        factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.axv[ip]/180.),2);
        factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.ayv[ip]/180.),2);
        factor += 1.;
        factor = TMath::Sqrt(factor);
        //      factor = 1.;///TEMP
        if(t_track.XGood(ip)){
          int iclx = t_track.GetClusterX_ID(ip);
          TrkCluster *clusterx = _trk_l1->GetCluster(iclx);
          float mip = TrkParams::GetMIP(clusterx->GetLadder()-1,clusterx->view-1);
          t_track.dedx_x[ip]   = clusterx->GetSignal()/(factor>0.?factor:1.)/(mip>0.?mip:1.); 
          t_track.multmaxx[ip] = clusterx->GetMultiplicity()*10000+clusterx->maxs;
          t_track.seedx[ip]    =  clusterx->clsignal[clusterx->indmax]/clusterx->clsigma[clusterx->indmax];//clusterx->GetSignalToNoise(0,7);
          t_track.xpu[ip]      = 0.;//complicato....
        }
        if(t_track.YGood(ip)){
          int icly = t_track.GetClusterY_ID(ip);
          TrkCluster *clustery = _trk_l1->GetCluster(icly);
          float mip = TrkParams::GetMIP(clustery->GetLadder()-1,clustery->view-1);
          t_track.dedx_y[ip]   = clustery->GetSignal()/(factor>0.?factor:1.)/(mip>0.?mip:1.); 
          t_track.multmaxy[ip] = clustery->GetMultiplicity()*10000+clustery->maxs;
          t_track.seedy[ip]    = clustery->clsignal[clustery->indmax]/clustery->clsigma[clustery->indmax];//clustery->GetSignalToNoise(0,7);
          t_track.ypu[ip]      = 0.;//complicato....
        }
      }
      ///////////////////////////////////////////////////////// end additional info

      //              cout << "ADD "<<endl;
      trackCandidates[mapIndex].insert(pair<int,TrkTrack>(t_track.GetNX()+t_track.GetNY(), TrkTrack(t_track)));
      if(t_track.GetNX()+t_track.GetNY()>5 && t_track.chi2<chi2Min) chi2Min=t_track.chi2;
      iCand++;
    }//end add-candidate condition
            
  }//end first loop on candidates
  if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;

  //=================================================================
  // ----------------------------------------------------------------
  // chi2 selection
  // ----------------------------------------------------------------
  //=================================================================
  pair <multimap<int,TrkTrack>::iterator, multimap<int,TrkTrack>::iterator> nxny[12];
  for(int inxny=0; inxny<12; inxny++){
    nxny[inxny] = trackCandidates[mapIndex].equal_range(inxny);
    for ( multimap<int,TrkTrack>::iterator it = nxny[inxny].first; it!=nxny[inxny].second; ++it){
      TrkTrack &trackC = (*it).second;
      //              cout << " NX+NY "<< trackC.GetNX()+trackC.GetNY() <<" chi2 "<<trackC.chi2<< " R(GV) "<<trackC.GetRigidity()<<endl;
      //            if( trackC.chi2==0 )continue;
      if( trackC.GetNX()+trackC.GetNY() > 5 && trackC.chi2>chi2Min )continue;
      trackCandidates[mapIndex+1].insert(pair<int,TrkTrack>(inxny, TrkTrack(trackC)));
    }
  }
  mapIndex++;

  if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
  
  //=================================================================
  // ----------------------------------------------------------------
  // save tracks
  // ----------------------------------------------------------------
  //=================================================================
  trkBegin = trackCandidates[mapIndex].begin();
  trkEnd   = trackCandidates[mapIndex].end();
  iCand = 0;
  TClonesArray &tracks = *_trkArray;
  tracks.Clear();
  for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){            
    TrkTrack trackCand = ((*trkIter).second);
    trackCand.seqno = iCand;
    //    trackCand.Dump();
    new(tracks[iCand]) TrkTrack(trackCand);
    iCand++;
  };


};
/**
 * Fill cluster map with tracker clusters
 */
void ExtTrkingAlg::FillClusterMap(multimap<int,int> &map,TrkLevel1* l1,Int_t vOffset){

  if(_debug)cout<<"-------- TRK CLUSTERs --------"<<endl;

  // _trk_cl.clear();

  if(!l1)return;

  ExtHit *trkCl = new ExtHit();

  for(int icl=0; icl<l1->nclstr(); icl++){

    // trkCl->Reset();

    // l1->GetCluster(icl)->GetLevel1Struct();// to use F77 routines //PERICOLOSO
    // trkCl->coordPU = (float)(l1->GetCluster(icl)->maxs) - 1. + l1->GetCluster(icl)->GetCOG();
    // trkCl->coordCm = 1000.; //TEMPORANEO forse non serve
    // trkCl->resCm = 1000.;//TEMPORANEO
    // trkCl->signal = l1->GetCluster(icl)->GetSignal();
    // trkCl->mult = l1->GetCluster(icl)->GetMultiplicity();
    // trkCl->view = l1->GetCluster(icl)->view-1;

    // _trk_cl.push_back(ExtHit(*trkCl));

    map.insert( pair<int,int>( l1->GetCluster(icl)->view-1 + vOffset, icl ));
    if(_debug)cout << icl << " - view(0-11) "<<l1->GetCluster(icl)->view-1<<" ladder(012) "<<l1->GetCluster(icl)->GetLadder()-1<<" strip "<<l1->GetCluster(icl)->maxs<<endl;
  }
  
  // l1->GetLevel1Struct(); //necessario per ritornare al contenuto originale del common 

  if(trkCl)delete trkCl;
};

/**
 * Fill cluster map with calorimeter clusters
 */
void ExtTrkingAlg::FillClusterMap(multimap<int,int> &map,CaloLevel1* l1,Int_t vOffset){

  if(_debug)cout<<"-------- CALO CLUSTERs -------- on VIEW < "<<_alg_nViewCal<<endl;

  if(_alg_nViewCal<=0)return;

  _cal_cl.clear();

  if(!l1)return;


  // ==========================
  // CALORIMETER CLUSTERING
  // ==========================
  CaloStrip st = CaloStrip(false);
  ExtHit *caloCl=NULL;
  Int_t nCl = 0;
  for(Int_t ih=0; ih<l1->istrip; ih++){//loop over calorimeter hits
    Int_t view = -1;  // 0-x 1-y
    Int_t plane = -1; // 0-21
    Int_t strip = -1; // 0-95
    Float_t mip = l1->DecodeEstrip(ih,view,plane,strip);
    
    // if(strip == 0 || strip == 95)cout <<" strip "<<strip<<" OK "<<endl;
    if(strip < 0 || strip > 95)cout <<" strip "<<strip<<" AHI AHI "<<endl;
    
    if(view<0)continue;
    if(plane<0)continue;
    if(strip<0)continue;
    if(mip<=0)continue;
    
    //      cout << ih << " -- "<<view<<" "<<plane<<" "<<strip<<endl;
    st.Set(view,plane,strip);// evaluate strip coordinates
    // 
    int   vv = plane*2+(int)(!view);           //0-43 odd x even y
    float cc = (view ? st.GetY() : st.GetX()); //coordinate
    float ss = mip;                            //signal
    float pp = (float)strip;                   //pitch units

    if( vv >= _alg_nViewCal ){//view not included
      if(!caloCl)continue; 
      if(caloCl->mult>0){ //if there is a cluster under construction...
        _cal_cl.push_back(ExtHit(*caloCl));//store it
        map.insert( pair<int,int>(caloCl->view + vOffset, nCl ));
        if(_debug)cout << nCl << " - view(0-43) "<< caloCl->view<<" strip "<<caloCl->start<<"-"<<caloCl->start+caloCl->mult-1<<endl;     
        nCl++;      
        caloCl->Reset();        
        caloCl->Set(strip,vv); // start new one... (without adding the strip)
      }
      continue;     //... and skip this view
    }    
    // --------------
    // @first hit:
    // --------------
    if(!caloCl){
      caloCl = new ExtHit(strip,vv);//create cluster starting 
      caloCl->Add(cc,pp,ss);//add strip
      continue;
    }
    //-------------------
    //check previous hit
    //-------------------     
    if(          
       vv != caloCl->view ||                  //hit on a different view
       strip > caloCl->start + caloCl->mult ||//not adjacent to previous one
       strip == 32 || //different ladder
       strip == 64 || //different ladder
       false ){       
      // is a new cluster:
      // store previous one, if not already stored
      if(caloCl->mult > 0){
        _cal_cl.push_back(ExtHit(*caloCl));
        map.insert( pair<int,int>(caloCl->view + vOffset, nCl ));
        if(_debug)cout << nCl << " - view(0-43) "<< caloCl->view<<" strip "<<caloCl->start<<"-"<<caloCl->start+caloCl->mult-1<<endl;     
        nCl++;      
      };
      //
      caloCl->Reset();      // reset
      caloCl->Set(strip,vv); // start new one from strip
      
    }
    caloCl->Add(cc,pp,ss); //add strip     
  };//end loop over calorimeter hits

  if(caloCl)delete caloCl;

  return;

};
/**
 * Evaluate positions
 */
bool ExtTrkingAlg::EvaluateClusterPosition_Tracker( int iclx, int icly,int lad,  int sen, float *xmABar, float *ymABar, float *zmAB){

  TrkParams::Load(1); 
  if( !TrkParams::IsLoaded(1) ){
    cout << "Bool_t ExtTrkingAlg::EvaluateClusterPosition_Tracker ---ERROR--- m.field not loaded "<<endl;
    return false;
  }    
  TrkParams::Load(4); 
  if( !TrkParams::IsLoaded(4) ){
    cout << "Bool_t ExtTrkingAlg::EvaluateClusterPosition_Tracker ---ERROR--- p.f.a. par. not loaded "<<endl;
    return false;
  }
  TrkParams::Load(5); 
  if( !TrkParams::IsLoaded(5) ){
    cout << "Bool_t ExtTrkingAlg::EvaluateClusterPosition_Tracker ---ERROR--- alignment par. not loaded "<<endl;
    return false;
  }
  //    cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;;
  int icx = iclx+1;//0=no-cluster,1-N
  int icy = icly+1;//0=no-cluster,1-N
  int sensor = sen+1;//<< convenzione "Paolo"
  int ip = (_trk_l1->GetCluster(TMath::Max(iclx,icly))->view-1)/2;//ip=0-5
  if(ip==5 && sensor!=0)sensor=3-sensor;//<< convenzione "Elena"
  int ladder = lad+1;
  float ax = (float)xmABar[3];
  float ay = (float)ymABar[3];
  float v[3];
  v[0]=xmABar[0];
  v[1]=ymABar[0];
  v[2]=zmAB[0];
  float bfx = 10*TrkParams::GetBX(v);//Tesla
  float bfy = 10*TrkParams::GetBY(v);//Tesla
  int ipp=ip+1;

  //  cout<<"xyzpam_("<<ipp<<","<<icx<<","<<icy<<","<<ladder<<","<<sensor<<"...)"<<endl;
  // e` una cosa bruttissima lo so ma che ci posso fare? mica posso fare tutto io!! 
  // ma lo sai a che ora me so` svegliata stamattina??
  xyzpam_(&ipp,&icx,&icy,&ladder,&sensor,&ax,&ay,&bfx,&bfy);
  //  cout<<" ip l s------------------- "<<ip<<ladder<<sensor<<" zm "<<track_.zm[ip]<<endl;
  extern cMini2track track_;
  
  xmABar[0]=track_.xm[ip];
  ymABar[0]=track_.ym[ip];
  zmAB[0]=track_.zm[ip];
  
  xmABar[1]=track_.xm_a[ip];
  ymABar[1]=track_.ym_a[ip];
  zmAB[1]=track_.zm_a[ip];//

  
  xmABar[2]=track_.xm_b[ip];
  ymABar[2]=track_.ym_b[ip];
  zmAB[2]=track_.zm_b[ip];//

  xmABar[4]=track_.resx[ip];
  ymABar[4]=track_.resy[ip];

  return true;

}
// bool ExtTrkingAlg::EvaluateClusterPosition_Calorimeter( int icl,  int sensor, float *xmABar, float *ymABar, float *zmAB, float def){

//   ExtHit cluster = _cal_cl[icl];
//   //  int plane = 6 + cluster.view;            //0...21
//   int plane = (int)(cluster.view / 2 );    //0...21
//   int ladder = (int)(cluster.coordPU / 32);//0...2
//   bool isX = (cluster.view)%2; 
//   bool isY = !isX;
  

//   bool ROTO = (sensor>=0);//if sensor is assigned (0,1,2) the strip can be rototraslated


//   CaloStrip st = CaloStrip(false);
//   //----------------------------
//   // Z-coordinate
//   //----------------------------
//   st.Set((int)isY,plane,ladder*32);//

//   float coordCm_Z0 = st.GetZ();
   
//   float size = 80.5; //(cm) larghezza del sensore + zone morte
//   //----------------------------
//   // strip-segment coordinates
//   //----------------------------
//   float coordCm_LA = cluster.coordCm;
//   float coordCm_LB = cluster.coordCm;
//   float coordCm_SA = (-1)*size*1.5; //full calorimeter size
//   float coordCm_SB = (+1)*size*1.5; //full calorimeter size
//   float coordCm_ZA = coordCm_Z0; 
//   float coordCm_ZB = coordCm_Z0; 

//   //----------------------------
//   // rototraslation
//   //----------------------------
//   if(ROTO){
//     //-----------------------------
//     // sensor origin coordinates
//     //-----------------------------
//     st.Set(   (int)isY, plane, ladder*32);//
//     float coordCm_L0 = ( isX ? st.GetX() : st.GetY());  
//     st.Set((int)(!isY), plane, sensor*32 );
//     float coordCm_S0 = (!isX ? st.GetX() : st.GetY());
    
//     coordCm_SA = (-1)*size*0.5 + (sensor-1)*size; //sensor size
//     coordCm_SB = (+1)*size*0.5 + (sensor-1)*size; //sensor size
    
//     //-----------------------------
//     // alignment pamameters
//     //-----------------------------
//     float fPitch = 1.;
//     float shift[]= {0.,0.,0.};//(L,S,Z)
//     float alpha  = 0.;
//     //-----------------------------

//     //coordinate relative al sensore
//     coordCm_LA = coordCm_LA - coordCm_L0; 
//     coordCm_SA = coordCm_SA - coordCm_S0;
//     coordCm_ZA = coordCm_ZA - coordCm_Z0;
//     coordCm_LB = coordCm_LB - coordCm_L0; 
//     coordCm_SB = coordCm_SB - coordCm_S0;
//     coordCm_ZB = coordCm_ZB - coordCm_Z0;

//     //variazione di pitch
//     coordCm_LA = coordCm_LA * fPitch;
//     coordCm_LB = coordCm_LB * fPitch;

//     //rotazione
//     float LA = coordCm_LA - alpha * coordCm_SA;
//     float SA = coordCm_SA + alpha * coordCm_LA;
//     float LB = coordCm_LB - alpha * coordCm_SB;
//     float SB = coordCm_SB + alpha * coordCm_LB;

//     //traslazione
//     coordCm_LA = LA         + shift[0];
//     coordCm_SA = SA         + shift[1];
//     coordCm_ZA = coordCm_ZA + shift[2];
//     coordCm_LB = LB         + shift[0];
//     coordCm_SB = SB         + shift[1];
//     coordCm_ZB = coordCm_ZB + shift[2];

//     //sistema di riferimento generale
    
//     coordCm_LA = coordCm_LA + coordCm_L0; 
//     coordCm_SA = coordCm_SA + coordCm_S0;
//     coordCm_ZA = coordCm_ZA + coordCm_Z0;
//     coordCm_LB = coordCm_LB + coordCm_L0; 
//     coordCm_SB = coordCm_SB + coordCm_S0;
//     coordCm_ZB = coordCm_ZB + coordCm_Z0;

//   }

//   xmABar[0] = -100.;
//   ymABar[0] = -100.;
//   zmAB[0]   = 0.5*coordCm_ZA+0.5*coordCm_ZB;
  
//   xmABar[1] = (isX ? coordCm_LA : coordCm_SA);
//   ymABar[1] = (isY ? coordCm_LA : coordCm_SA);
//   zmAB[1]   = coordCm_ZA;

  
//   xmABar[2]= (isX ? coordCm_LB : coordCm_SB);
//   ymABar[2]= (isY ? coordCm_LB : coordCm_SB);
//   zmAB[2]  = coordCm_ZB;


//   //----------------------------
//   // evaluate spatial resolution
//   //----------------------------

 
//   float res_dig = 0.24/sqrt(12.);

//   float res_ms = 0.;
  
//   if(def!=0.){
    
//     float rig = 1/def; if(rig<0)rig = -rig;
//     float factor = 0.;
//     factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*xmABar[3]/180.),2);
//     factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*ymABar[3]/180.),2);
//     factor += 1.;
//     factor = TMath::Sqrt(factor);
    
//     int view = cluster.view; //0... 43

//     int nW = (int)((view + 1)/2);
//     float dW = 0.74*factor;//X0
//     float dW_cm = 0.26*factor;//cm
//     float dSi_cm = 0.; //cm
//     if(isY && plane>0) dSi_cm = ( plane%2 ? 0.228 : 0.428);
//     dSi_cm *= factor;
//     // multiple scattering angle
//     float theta = 0.;//ms
//     theta = 0.0136/rig;
//     theta = theta * sqrt(dW) * (1.+0.038*log(dW));
//     float Vt = theta*theta;            //cov(theta,theta)
//     float Vy = Vt * dW_cm * dW_cm /3.; //cov(y,y)
//     float Vty = 0.87 * Vt * Vy;        //cov(theta,y)
//     float Vc = Vy + dSi_cm * dSi_cm * Vt + 2. * dSi_cm * Vty;
//     float rms = sqrt( Vc );    
//     // if(isY)t_track.resy[ip] = sqrt(rms*rms + t_track.resy[ip]*t_track.resy[ip]);
//     // else   t_track.resx[ip] = sqrt(rms*rms + t_track.resx[ip]*t_track.resx[ip]);
//     // if(isY)          VMS[nW] = rms*rms;
//     // if(_debug)cout <<endl<<view<<" nW "<<nW<<" res(MS) "<<rms<<" factor "<<factor;
      
//     if(view>=3)cout << " --> NB error computation not accurate ";
//     res_ms = rms;

//   }


//   xmABar[4]= (isX ? sqrt(res_dig*res_dig+res_ms*res_ms) : 1000.);
//   ymABar[4]= (isY ? sqrt(res_dig*res_dig+res_ms*res_ms) : 1000.);


// }
bool ExtTrkingAlg::EvaluateClusterPosition_Calorimeter( int icl, int& sisensor, float *xmABar, float *ymABar, float *zmAB, float def){

  ExtHit cluster = _cal_cl[icl];
  //  int plane = 6 + cluster.view;            //0...21
  int plane = (int)(cluster.view / 2 );    //0...21
  int ladder = (int)(cluster.coordPU / 32);//0...2
  //  int view = cluster.view;//0...43    // EM, unused variable
  bool isX = (cluster.view)%2; 
  bool isY = !isX;
  

//  int sisensor = -1;

  
//  CaloStripRoto temp;
  
  bool ROTO = ( xmABar[0]!=0. && ymABar[0]!=0.);
// ------------------------------------------------------------------------------
//
// check if the cluster belongs to the hit sensor (xmABar[0],ymABar[0] must be provided)
// evaluate the hit sensor number
//
// ------------------------------------------------------------------------------
  if(ROTO){


      for(int is=0;is<3; is++){
          int sis = GetCaloSiSensor((int)isY,ladder,is);
//        cout << "view, ladder, is :"<<(int)isY<<ladder<<is<<" ---> sis "<<sis<<endl;
          int index = (plane*2+1-(int)isY)*9 + sis;
          CaloStripRoto extst = _caloStripRoto[index];
//        cout << "index "<<index<<endl;
//        cout << " -> view     "<<extst.GetView()<<endl;
//        cout << " -> plane    "<<extst.GetPlane()<<endl;
//        cout << " -> sisensor "<<extst.GetSiSensor()<<endl;
//        cout << " -> ladder   "<<extst.GetLadder(<<" "<<GetCaloLadder(extst.GetView(),extst.GetSiSensor())<<endl;       
          bool inside = extst.SensorContains(xmABar[0],ymABar[0]);
//        cout << " inside "<<sis<<" ?? "<<inside<<endl;
          if( inside ){
              sisensor = sis;
              continue;
          }
      }

      if(sisensor==-1){
//        cout << "ATTENZIONE!!! sisensor "<<sisensor<<" x,y "<<xmABar[0]<<" "<<ymABar[0]<<endl;
      }
  }

// ------------------------------------------------------------------------------
//
// if no sensor is hit, skip rototraslation
//
// ------------------------------------------------------------------------------
  if(sisensor==-1)ROTO=false;
// ------------------------------------------------------------------------------
//
// now evaluate strip A-B coordinates
//
// ------------------------------------------------------------------------------
  if(ROTO){

      int index = (plane*2+1-(int)isY)*9 + sisensor;
      CaloStripRoto extst = _caloStripRoto[index];
//      CaloStrip st = extst.st;
      
//       cout << " extst.Set("<<(int)isY<<", "<<plane<<", "<<cluster.coordPU<<", "<<sensor<<") "<<endl;
//       extst.Set((int)isY, plane, cluster.coordPU, sensor);

      if( sisensor != extst.GetSiSensor() )cout << " extst.GetSiSensor() "<<extst.GetSiSensor()<<" != "<<sisensor<<endl;
      if( plane != extst.GetPlane() )cout << " extst.GetPlane() "<<extst.GetPlane()<<" != "<<plane<<endl;
      if( (int)isY != extst.GetView() )cout << " extst.GetView() "<<extst.GetView()<<" != "<<isY<<endl;
//      cout << " extst.SetStrip("<<cluster.coordPU<<") "<<endl;
      extst.SetStrip(cluster.coordPU);
      

      xmABar[0] = -100.;
      ymABar[0] = -100.;
      zmAB[0]   = 0.5*extst.fZA+0.5*extst.fZB;
      
      xmABar[1] = extst.fXA;
      ymABar[1] = extst.fYA;
      zmAB[1]   = extst.fZA;
      
      
      xmABar[2]= extst.fXB;
      ymABar[2]= extst.fYB;
      zmAB[2]  = extst.fZB;

      xmABar[4]= (isX ? extst.GetSpatialResolution(def, xmABar[3], ymABar[3],  1.) : 1000.);
      ymABar[4]= (isY ? extst.GetSpatialResolution(def, xmABar[3], ymABar[3],  1.) : 1000.);
  

  }else{
      

      CaloStrip st = CaloStrip(false);
      //----------------------------
      // Z-coordinate
      //----------------------------
      st.Set((int)isY,plane,ladder*32);//

      float coordCm_Z0 = st.GetZ();
   
      float size = 80.5; //(cm) larghezza del sensore + zone morte
      //----------------------------
      // strip-segment coordinates
      //----------------------------
      float coordCm_LA = cluster.coordCm;
      float coordCm_LB = cluster.coordCm;
      float coordCm_SA = (-1)*size*1.5; //full calorimeter size
      float coordCm_SB = (+1)*size*1.5; //full calorimeter size
      float coordCm_ZA = coordCm_Z0; 
      float coordCm_ZB = coordCm_Z0; 

      //----------------------------
      // rototraslation
      //----------------------------

      xmABar[0] = -100.;
      ymABar[0] = -100.;
      zmAB[0]   = 0.5*coordCm_ZA+0.5*coordCm_ZB;
  
      xmABar[1] = (isX ? coordCm_LA : coordCm_SA);
      ymABar[1] = (isY ? coordCm_LA : coordCm_SA);
      zmAB[1]   = coordCm_ZA;

  
      xmABar[2]= (isX ? coordCm_LB : coordCm_SB);
      ymABar[2]= (isY ? coordCm_LB : coordCm_SB);
      zmAB[2]  = coordCm_ZB;

  
      //----------------------------
      // evaluate spatial resolution
      //----------------------------

 
      float res_dig = 0.24/sqrt(12.);

      float res_ms = 0.;
  
      if(def!=0.){
    
          float rig = 1/def; if(rig<0)rig = -rig;
          float factor = 0.;
          factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*xmABar[3]/180.),2);
          factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*ymABar[3]/180.),2);
          factor += 1.;
          factor = TMath::Sqrt(factor);
    
          int view = cluster.view; //0... 43

          //      int nW = (int)((view + 1)/2); // EM, unused variable
          float dW = 0.74*factor;//X0
          float dW_cm = 0.26*factor;//cm
          float dSi_cm = 0.; //cm
          if(isY && plane>0) dSi_cm = ( plane%2 ? 0.228 : 0.428);
          dSi_cm *= factor;
          // multiple scattering angle
          float theta = 0.;//ms
          theta = 0.0136/rig;
          theta = theta * sqrt(dW) * (1.+0.038*log(dW));
          float Vt = theta*theta;            //cov(theta,theta)
          float Vy = Vt * dW_cm * dW_cm /3.; //cov(y,y)
          float Vty = 0.87 * Vt * Vy;        //cov(theta,y)
          float Vc = Vy + dSi_cm * dSi_cm * Vt + 2. * dSi_cm * Vty;
          float rms = sqrt( Vc );         
          // if(isY)t_track.resy[ip] = sqrt(rms*rms + t_track.resy[ip]*t_track.resy[ip]);
          // else   t_track.resx[ip] = sqrt(rms*rms + t_track.resx[ip]*t_track.resx[ip]);
          // if(isY)          VMS[nW] = rms*rms;
          // if(_debug)cout <<endl<<view<<" nW "<<nW<<" res(MS) "<<rms<<" factor "<<factor;
      
          if(view>=3)cout << " --> NB error computation not accurate ";
          res_ms = rms;

      }

      xmABar[4]= (isX ? sqrt(res_dig*res_dig+res_ms*res_ms) : 1000.);
      ymABar[4]= (isY ? sqrt(res_dig*res_dig+res_ms*res_ms) : 1000.);

  }

  return true; // EM, warning: no return statement in function returning non-void [-Wreturn-type]
}


/**
 * Process the event
 */
void ExtTrkingAlg::ProcessEvent1(Bool_t force){

  if(_debug)cout << " |---------------------------------------------------| "<<endl;
  if(_debug)cout << " void ExtTrkingAlg::ProcessEvent1("<<force<<") "<<endl;

  
  if(_debug && !_trk_l1)cout << " Missing TrkLevel1 object --- ciao! "<<endl;
  if(!_trk_l1)return;


  //===========================
  // pre-selection  
  //===========================
  
  Bool_t RETRACK = CheckEvent();

  if(_debug)cout <<" Pre-selection result: "<<RETRACK<<endl;
  if( !force && !RETRACK )return;
  

  //===========================
  // tracking
  //===========================

  // TClonesArray &tracks = *_trkArray;
  // tracks.Clear();
  TClonesArray &tracks = *_trkArray;
  tracks.Clear();


  multimap<int,ExtTrack> trackCandidates[NEXTVIEWS]; // map: last hit view vs TrkTrack
  int mapIndex = 0;
  multimap<int,ExtTrack>::iterator trkBegin;
  multimap<int,ExtTrack>::iterator trkEnd;
  int iCand;

  // ------------------------------------------------
  // fill a map VIEW-CLUSTER
  // ------------------------------------------------
  multimap<int,int> clusterMap;
  FillClusterMap(clusterMap,_trk_l1,0);// trk clusters    
  FillClusterMap(clusterMap,_cal_l1,12);// calorimeter clusters
  

  // ------------------------------------------------
  // define iterators over clusters on the same view
  // ------------------------------------------------
  pair <multimap<int,int>::iterator, multimap<int,int>::iterator> ret[NEXTVIEWS];
  multimap<int,int>::iterator it[NEXTVIEWS];
  for(int iv=0; iv<NEXTVIEWS; iv++)ret[iv] = clusterMap.equal_range(iv);

  
  //====================================================================
  // PARAMETERS
  //====================================================================
  // ------------------------------------------------
  // number of points per candidate
  // ------------------------------------------------
  int nClFixX =  _alg_nClFixX; // n.hits required on X view
  int nClFixY = _alg_nClFixY; // n.hits required on X view
  int nTrackMAX =  _alg_nTrackMAX; // 
  //==================================================
  // ------------------------------------------------
  // start: one candidate for each cluster
  // ------------------------------------------------
  //==================================================
  if(_debug)cout<<"-------- MINIMAL TRACK-CANDIDATES --------";
  if(_debug)cout<<nClFixX<<"X"<<nClFixY<<"Y"<<endl;
  for(int iv=0; iv<NEXTVIEWS; iv++){                                      //loop over all views
    
    if( !(iv%2) && (int)(0.5*iv) >  NEXTVIEWS/2-nClFixX)continue;
    if(  (iv%2) && (int)(0.5*iv) >  NEXTVIEWS/2-nClFixY)continue;

    if(iv >= NEXTVIEWS - _alg_nViewCal )continue; //skip calo hits 
    
    for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view


      int icl = it[iv]->second; 

      TrkCluster *cluster = _trk_l1->GetCluster(icl);
      //      cluster->Dump();
      int ladder = cluster->GetLadder()-1; //
      int plane = cluster->GetPlane()-1;
      bool isY = (cluster->view)%2;// F77 1-12 ??
      bool isX = !isY;
      

      /// TRACKER ==============================================================================
      for (int is=0; is<2; is++){                                  //loop over sensors

        //cout <<endl<< " view "<<iv<<" plane "<<plane<<" cluster "<<icl<<" on sensor "<<is<<endl;
        
        ExtTrack& t_track =  *_extTrack;

        t_track.ResetXY();
        for(int ip =0 ; ip<t_track.nplanes; ip++ )t_track.SetZ(ip,_zMech[ip]);
        //      t_track.Dump();


        // >>>> insert point
        float xmABar[]={0.,0.,0.,0.,0.}; // inizializzare a zero l'angolo, come prima stima per il PFA
        float ymABar[]={0.,0.,0.,0.,0.};
        float zmAB[]={t_track.zm[plane],0.,0.};//mechanical z-coordinate
        EvaluateClusterPosition_Tracker( (isX?icl:-1), (isY?icl:-1), ladder, is, xmABar, ymABar, zmAB); 
        
        if( isX ){
          t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
          t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
          t_track.SetXGood(plane,icl+1,ladder,is);      
          //      t_track.SetYGood(plane,0,ladder,is);  
        }
        if( isY ){
          t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
          t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
          t_track.SetYGood(plane,icl+1,ladder,is);
          //      t_track.SetXGood(plane,0,ladder,is);
        }
        trackCandidates[mapIndex].insert(pair<int,ExtTrack>(iv, t_track));

      }
      

    }
  }

  //==============================================================
  // -------------------------------------------------------------
  // next: add views to each candidate, up to nClFixX+nClFixY
  // -------------------------------------------------------------
  //==============================================================
  for( ;mapIndex < nClFixX+nClFixY-1 ;++mapIndex){

    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    if(_debug)cout<<"MAP "<<mapIndex<<" ----------- n.candidates "<<trackCandidates[mapIndex].size()<<endl;

    for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){


      ExtTrack trackCand = (*trkIter).second;
      int lastView = (*trkIter).first; //last inserted view
      bool lastIsTrk = (lastView < NEXTVIEWS - _alg_nViewCal);
      int lastPlane = (int)(0.5*lastView);
      if(!lastIsTrk)lastPlane =lastView-6;
      int lastLadder = trackCand.GetLadder(lastPlane);
      int lastSensor = trackCand.GetSensor(lastPlane);

      
      // if(_debug)cout<<"MAP "<<mapIndex<<" last v "<<lastView<<" NX "<<trackCand.GetNX()<<" NY "<<trackCand.GetNY()<<endl;

      for(int iv=lastView+1; iv<NEXTVIEWS; iv++){                             //loop over next views

        bool isTrk = (iv < NEXTVIEWS - _alg_nViewCal);
        bool isX = (iv)%2;
        bool isY = !isX;
        
                     
        // if(  (iv%2) && (int)(0.5*iv) > 6-nClFixX+trackCand.GetNX())continue; //not enough x-view left
        // if( !(iv%2) && (int)(0.5*iv) > 6-nClFixY+trackCand.GetNY())continue; //not enough y-view left

        // if(  (iv%2) && trackCand.GetNX()==nClFixX )continue;//ok?
        // if( !(iv%2) && trackCand.GetNY()==nClFixY )continue;//ok?

        if( isX && trackCand.GetNX()==nClFixX )continue;//ok?
        if( isY && trackCand.GetNY()==nClFixY )continue;//ok?


        if( isX && (int)(0.5*iv) >  6+_alg_nViewCal/2     - nClFixX+trackCand.GetNX())continue; //not enough x-view left
        if( isY && (int)(0.5*iv) >  6+(_alg_nViewCal+1)/2 - nClFixY+trackCand.GetNY())continue; //not enough y-view left


        for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view

          int icl = it[iv]->second;

          if(isTrk){
            /// TRACKER ==================================================
            

            TrkCluster *cluster = _trk_l1->GetCluster(icl);
            int ladder = cluster->GetLadder()-1;//(_trk_l1->GetCluster(icl)->maxs-1)%1024 + 1;
            int plane = cluster->GetPlane()-1;
            // bool isY = (cluster->view)%2;
            // bool isX = !isY;
            
            if( plane==lastPlane && ladder!=lastLadder)continue;
            

            for (int is=0; is<2; is++){                                  //loop over sensors
                  
              if( plane==lastPlane && is!=lastSensor)continue;
              
              // if(_debug)cout<<"MAP "<<mapIndex<<" last v "<<lastView;
              // if(_debug)cout<<" +cl "<<icl<<" s"<<is<<endl;
              
              ExtTrack t_track = ExtTrack();
              trackCand.Copy(t_track);                             //copy the candidate parameters...
              
              bool isXY = (isX && t_track.YGood(plane)) || (isY && t_track.XGood(plane)) ;
              

              float xmABar[]={0.,0.,0.,0.,0.}; //importante inizializzare a zero, come prima stima per il PFA
              float ymABar[]={0.,0.,0.,0.,0.};
              float zmAB[]={t_track.zm[plane],0.,0.};

              if( isX ){
                int iclx = icl;
                int icly = (isXY ? t_track.GetClusterY_ID(plane) : -1 );

                EvaluateClusterPosition_Tracker( iclx, icly , ladder, is, xmABar, ymABar, zmAB); 

                if(isXY){
                  t_track.SetXY(plane,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
                  t_track.SetZ(plane,zmAB[0]);
                }else{
                  t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
                  t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
                }
                t_track.SetXGood(plane,iclx+1,ladder,is);         //...add a point...
                t_track.SetYGood(plane,icly+1,ladder,is);
              }
              //
              if( isY ){
                int iclx = (isXY ? t_track.GetClusterX_ID(plane) : -1 );
                int icly = icl;

                EvaluateClusterPosition_Tracker( iclx, icly , ladder, is, xmABar, ymABar, zmAB); 

                if(isXY){
                  t_track.SetXY(plane,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
                  t_track.SetZ(plane,zmAB[0]);
                }else{
                  t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
                  t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
                };
                t_track.SetYGood(plane,icly+1,ladder,is);
                t_track.SetXGood(plane,iclx+1,ladder,is);
              }

              trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(iv, t_track)); //...and store a new candidate
            };//end loop over sensors
            
          }else{              
            /// CALORIMETER =================================================
            ExtHit cluster = _cal_cl[icl];
            int plane = 6 + cluster.view;

            if( plane==lastPlane )continue;

            int ladder = (int)(cluster.coordPU / 32);
            int sensor = -1; //not yet known
            bool isX = (cluster.view)%2;
            bool isY = !isX;

            ExtTrack t_track = ExtTrack();
            trackCand.Copy(t_track);   

            float xmABar[]={0.,0.,0.,0.,0.}; //importante inizializzare a zero, come prima stima per il PFA
            float ymABar[]={0.,0.,0.,0.,0.}; // sensor is not evaluated
            float zmAB[]={0.,0.,0.};
            EvaluateClusterPosition_Calorimeter( icl, sensor, xmABar, ymABar, zmAB); 


            if( isX ){
              t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
              t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
              t_track.SetXGood(plane,icl+1, ladder, sensor);         //...add a point.
              t_track.SetYGood(plane,    0, ladder, sensor);     //...add a point...
            }      
            if( isY ){
              t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
              t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
              t_track.SetYGood(plane,icl+1, ladder, sensor);
              t_track.SetXGood(plane,    0, ladder, sensor);
            }
            trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(iv, t_track)); //...and store a new ca

          }

          if( trackCandidates[mapIndex].size() > (uint)nTrackMAX ){ // EM, compilation warning comparison between signed and unsigned
            if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<" > "<<nTrackMAX<<endl;
            return;//to many candidates
          }

        };//end loop over clusters on the view iv
      };//endl loop over views
    }//end loop over candidates
  }//end iterations
  

  // if( trackCandidates[mapIndex].size() > nTrakMAX ){
  //   if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<" > "<<nTrackMAX<<endl;
  //   return;//to many candidates
  // }

  // if(_debug){
  //   trkBegin = trackCandidates[mapIndex].begin();
  //   trkEnd   = trackCandidates[mapIndex].end();
  //   for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
  //     ExtTrack &trackCand = (*trkIter).second;
  //     cout << " >> ";
  //     cout << " X: ";
  //     for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
  //     cout << " Y: ";
  //     for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
        
  //     cout << endl;
  //   }
  // }

  if(_debug)cout<<"MAP "<<mapIndex<<" ----------- n.candidates "<<trackCandidates[mapIndex].size()<<endl;

  if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;

  if(trackCandidates[mapIndex].size()==0)return;//no good candidates

  //===============================================================================================
  // -------------------------------------------------------------
  // first selection of track candidates, based on chi2
  // -------------------------------------------------------------
  //===============================================================================================
  if(_debug)cout<<"-------- CHI2 --------"<<endl;
  trkBegin = trackCandidates[mapIndex].begin();
  trkEnd   = trackCandidates[mapIndex].end();
  mapIndex++;
  iCand = 0;
  for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
    int lastView = (*trkIter).first;
    ExtTrack &trackCand = (*trkIter).second;      
    int ifail=0;
    trackCand.ResetFit();
    trackCand.Fit(0.,ifail,0);
    if(ifail!=0)trackCand.ResetFit();

    // -----------------------------------
    // first selection of track candidates
    // -----------------------------------
    if(TMath::IsNaN(trackCand.chi2))continue;
    if(trackCand.chi2 <= 0. || trackCand.chi2 > 1.e7)continue; 
    //    if(trackCand.chi2 < 0. || trackCand.chi2 > 10. )continue; 
    iCand++;    
          
    if(_debug){     
     cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
      cout <<" | al: ";
      for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
      cout <<" | ";
      cout << " X: ";
      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
      // for(int ip=0; ip<trackCand.nplanes; ip++)
      // if(trackCand.GetClusterX_ID(ip)>=0)cout << trackCand.GetClusterX_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";
      cout << " | ";
      cout << " Y: ";
      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
      cout << endl;
      // for(int ip=0; ip<trackCand.nplanes; ip++)
      //        if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";       cout << endl;
    }

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // evaluated coordinates (to define GF)
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // int    ngf = TrkParams::nGF;
    // float *zgf = TrkParams::zGF;
    // Trajectory tgf = Trajectory(ngf,zgf);
    // _tgf->DoTrack(trackCand.al);//<<<< integrate the trajectory
    // for(int ip=0; ip<ngf; ip++){
    //   trackCand.xGF[ip] = tgf->x[ip];
    //   trackCand.yGF[ip] = tgf->y[ip];
    // }
        
    trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, trackCand));

  }
  if(trackCandidates[mapIndex].size()==0)return;//no good candidates

  if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
  //=================================================================
  // ----------------------------------------------------------------
  // CALORIMETER
  // ----------------------------------------------------------------
  //=================================================================
  if(_cal_l1){
    if(_debug)cout<<"-------- CALORIMETER --------"<<endl;
    // -----------------------------------
    // evaluate calorimeter variables
    // -----------------------------------
    // float caloCoord[2][22][96];
    // float caloZ[44];
    // for(int view=0; view<2; view++){
    //   for(int plane=0; plane<22; plane++){
    //  for(int strip=0; strip<96; strip++){
    //    CaloStrip st = CaloStrip(false);
    //    st.Set(view,plane,strip);
    //    caloCoord[view][plane][strip]=(view ? st.GetY() : st.GetX());
    //    caloZ[plane*2+(view?0:1)]=st.GetZ();
    //  }
    //   }
    // }
    // Trajectory caloTj = Trajectory(44,caloZ);
    //  for(int iz=0; iz<44;iz++)cout<<caloZ[iz]<<endl;
    //test calo
    //    float caloEvent[2][22][96];
    CaloLevel1* l1=_cal_l1;//event->GetCaloLevel1();
    //    if(!l1)return;
    vector<float> caloChi2;
    vector<int> caloHit; int caloHitMAX=0;
    vector<float> caloMip;
    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){        
      ExtTrack trackCand = ((*trkIter).second);   
            
      // trackCand.DoTrack(&caloTj);//integrate the trajectory through the calorimeter
      _caloTj->DoTrack(trackCand.al,trackCand.zini);
      //      _caloTj->Dump();
      
      float chi2 = 0.;
      int  nhit = 0;
      float cmip = 0.;
      for(Int_t ih=0; ih<l1->istrip; ih++){
        Int_t view = -1;
        Int_t plane = -1;
        Int_t strip = -1;
        Float_t mip = l1->DecodeEstrip(ih,view,plane,strip);
        if(view<0)continue;
        if(plane<0)continue;
        if(strip<0)continue;
        if(mip<=0)continue;
        float dxy = _caloCoord[view][plane][strip]-(view?_caloTj->y[plane*2+(view?0:1)]:_caloTj->x[plane*2+(view?0:1)]);
        if(plane<11) chi2+= mip * TMath::Power(dxy,2.);
        //            caloZ[plane*2+(view?0:1)]=st.GetZ();
        if( fabs(dxy) < 2. && plane<11)nhit++;
        if( fabs(dxy) < 2. && plane<11)cmip+=mip;
      };
      if(l1->istrip>0)chi2 = chi2/l1->istrip;
            
      caloHit.push_back(nhit); 
      if(nhit>caloHitMAX)caloHitMAX=nhit;

      if(_debug){           
        cout << " CALO ";
        //      cout << " chi2 "<< chi2 <<" ";
        cout << " nhit (< 11th plane) "<< nhit <<" ";
        //      cout << " mip  "<< cmip <<" ";
        cout <<endl;
      };
    }

          
    //=================================================================
    // ----------------------------------------------------------------
    // selection of track candidates, based on n.calorimeter hits
    // ----------------------------------------------------------------
    //=================================================================
    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    mapIndex++;
    iCand = 0;
    for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
      int lastView = (*trkIter).first;
      ExtTrack &trackCand = (*trkIter).second;
      if(caloHit[iCand]==caloHitMAX)trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, trackCand));
      iCand++;
    }
    if(trackCandidates[mapIndex].size()==0)return;//no good candidates
    if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;

  }


  //=================================================================
  // ----------------------------------------------------------------
  // TOF
  // ----------------------------------------------------------------
  //=================================================================
  if(_tof_l2){
    if(_debug)cout<<"-------- TOF --------"<<endl;
    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    ToFLevel2* tl2=_tof_l2;//event->GetToFLevel2(); 
    //   if(!tl2)return;
    iCand =0;
    vector<int> tofHit; int tofHitMAX=0; 
    for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
      //            int lastView = (*trkIter).first;
      ExtTrack &trackCand = (*trkIter).second;

      _tgf->DoTrack(trackCand.al);//<<<< integrate the trajectory

       for(int ip=0; ip<TrkParams::nGF; ip++){
        trackCand.xGF[ip] = _tgf->x[ip];
        trackCand.yGF[ip] = _tgf->y[ip];
      }

      int nhit =0;
      for(int iToF=0; iToF<6; iToF++){
              
        int iGF = iToF;
        if(iToF>3)iGF = iToF+8;
        int iPaddle = tl2->GetPaddleIdOfTrack( trackCand.xGF[iGF], trackCand.yGF[iGF], iToF);
        if( tl2->HitPaddle(iToF,iPaddle) )nhit++;
              
      }
      iCand++;
            
      tofHit.push_back(nhit); 
      if(nhit>tofHitMAX)tofHitMAX=nhit;

      if(_debug){
        cout << " TOF nhit "<< nhit <<" ";
        cout <<endl;
      }
            
    }
    //=================================================================
    // ----------------------------------------------------------------
    // selection of track candidates, based on n.tof hits
    // ----------------------------------------------------------------
    //=================================================================
    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    mapIndex++;
    iCand = 0;
    for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
      int lastView = (*trkIter).first;
      ExtTrack &trackCand = (*trkIter).second;
      if(tofHit[iCand]==tofHitMAX)trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, ExtTrack(trackCand)));
      iCand++;

    }
    if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
    if(trackCandidates[mapIndex].size()==0)return;//no good candidates

  }

  if(_debug){     
    
    cout<< " Minimal-track SUMMARY:"<<endl;

    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
      ExtTrack &trackCand = (*trkIter).second;
      cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
      cout <<" | al: ";
      for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
      cout <<" | ";
      cout << " X: ";
      //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
      for(int ip=0; ip<trackCand.nplanes; ip++){
        if(trackCand.GetClusterX_ID(ip)>=0)
          cout << setw(2)<<trackCand.GetClusterX_ID(ip);
        else
          cout << "__";
      }
      cout << " | ";
      cout << " Y: ";
      //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
      for(int ip=0; ip<trackCand.nplanes; ip++){
        if(trackCand.GetClusterY_ID(ip)>=0)
          cout << setw(2)<<trackCand.GetClusterY_ID(ip);
        else
          cout << "__";
      }
      cout << " S:";for(int ip=0; ip<trackCand.nplanes; ip++)cout<< trackCand.GetSensor(ip)+1;
      cout << endl;

    }
  }
  
  //=================================================================
  // ----------------------------------------------------------------
  // TRACK MERGING
  // ----------------------------------------------------------------
  //=================================================================
  if(_debug)cout<<"-------- MERGING --------";
  trkBegin = trackCandidates[mapIndex].begin();
  trkEnd   = trackCandidates[mapIndex].end();
  mapIndex++;
  float chi2Min=1e10;
  iCand = 0;
  for( multimap<int,ExtTrack>::iterator trkIter1 = trkBegin; trkIter1!=trkEnd; ++trkIter1){ //loop over tracks (1)
    
    ExtTrack &trackCand1 = (*trkIter1).second;// get candidate 1...
            
    ExtTrack t_track = ExtTrack();// ...create a new ExtTrack instance...
    trackCand1.Copy(t_track);     // ...and make a copy
                   
    bool HBM = false; //has been merged

    for( multimap<int,ExtTrack>::iterator trkIter2 = trkBegin; trkIter2!=trkEnd; ++trkIter2){//loop over tracks (2)
              
      ExtTrack &trackCand2 = (*trkIter2).second; // get candidate 2...

      bool CBM = true;
      for(int ip=0; ip<NEXTPLANES; ip++){ //loop over planes
        if( t_track.XGood(ip)  &&  
            trackCand2.XGood(ip) &&  
            ( t_track.GetClusterX_ID(ip)!=trackCand2.GetClusterX_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
            true )CBM=false; // if both have a x-hit and it is not the same ---> they cannot be merged
        if(!CBM)break;       
        if( t_track.YGood(ip)  && 
            trackCand2.YGood(ip) && 
            ( t_track.GetClusterY_ID(ip)!=trackCand2.GetClusterY_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
            true )CBM=false; // if both have a y-hit and it is not the same ---> they cannot be merged
        if(!CBM)break;        
      }
      //      if(_debug)cout << " can be merged ? "<<CBM<<endl;
      if(!CBM)continue; //go to the next candidate
      ///////////
      //merging//
      ///////////    
      // cout <<endl<<" check merging al (1) ";
      // for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
      // cout <<endl<<"       with    al (2) ";
      // for(int i=0; i<5; i++)cout <<setw(10)<< trackCand2.al[i];
  
      for(int ip=0; ip<NEXTPLANES; ip++){             
        if( !(t_track.XGood(ip)) &&  trackCand2.XGood(ip) ){      
          int ygood_temp = t_track.ygood[ip];
          
          if(t_track.YGood(ip)){ 
            t_track.SetXY(ip,trackCand2.xm[ip],t_track.ym[ip],trackCand2.resx[ip],t_track.resy[ip]);
            //      cout<<"confronta "<<trackCand2.zm[ip]<<" - "<<t_track.zm[ip]<<" = "<<trackCand2.zm[ip]-t_track.zm[ip]<<endl;
          }else{ 
            t_track.SetX(ip,trackCand2.xma[ip],trackCand2.xmb[ip],trackCand2.yma[ip],trackCand2.ymb[ip],trackCand2.resx[ip]);
            t_track.SetZ(ip,0.5*(trackCand2.zma[ip]+trackCand2.zmb[ip]));
          }
          t_track.xgood[ip] = trackCand2.xgood[ip];//assign cluster + sensor
          t_track.ygood[ip] = ygood_temp;//assign cluster + sensor

          HBM = true;    
        }               
        if( !(t_track.YGood(ip)) &&  trackCand2.YGood(ip) ) {     

          int xgood_temp = t_track.xgood[ip];
          
          if(t_track.XGood(ip)){ 
            t_track.SetXY(ip,t_track.xm[ip],trackCand2.ym[ip],t_track.resx[ip],trackCand2.resy[ip]);     
            //      cout<<"confronta "<<trackCand2.zm[ip]<<" - "<<t_track.zm[ip]<<" = "<<trackCand2.zm[ip]-t_track.zm[ip]<<endl;
          }else{ 
            t_track.SetY(ip,trackCand2.xma[ip],trackCand2.xmb[ip],trackCand2.yma[ip],trackCand2.ymb[ip],trackCand2.resy[ip]);
            t_track.SetZ(ip,0.5*(trackCand2.zma[ip]+trackCand2.zmb[ip]));
          }
          t_track.ygood[ip] = trackCand2.ygood[ip];//assign cluster + sensor 
          t_track.xgood[ip] = xgood_temp;//assign cluster + sensor
          HBM = true;
        }
      };         
      
    }//end loop over tracks (2)


    // cout << "<<  ";
    // for(int ip=0; ip<6; ip++) cout <<setw(4)<< t_track.GetClusterX_ID(ip);
    // for(int ip=0; ip<6; ip++) cout <<setw(4)<< t_track.GetClusterY_ID(ip);
    // cout << endl;
    // cout<< " TRK chi2 "<<setw(10)<<t_track.chi2;
    // cout <<" | al ";
    // for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
    // cout << endl;
 
    ///////////////////////////////////////////////////////// FIT START
    // int ifail=0;
    // //           TrkParams::SetDebugMode();
    // t_track.FitReset();       //angles set to zero
    // t_track.Fit(0.,ifail,0,1);//re-evaluate cluster position from level1
    // if( TMath::IsNaN(t_track.chi2) )continue;
    // if( ifail != 0 )continue;
    // t_track.Fit(0.,ifail,0,1);//evaluate again taking into account the track angle
    // if( TMath::IsNaN(t_track.chi2) )continue;
    // if( ifail != 0 )continue;
    ///////////////////////////////////////////////////////// FIT FINISH
    // cout << "Merged: "<<endl;
    // cout << "X "; for(int ip=0; ip<6; ip++)cout << t_track.XGood(ip); cout << endl;
    // cout << "Y "; for(int ip=0; ip<6; ip++)cout << t_track.YGood(ip); cout << endl;
   //
    //check if track has been already inserted
    //
    multimap<int,ExtTrack>::iterator trkB = trackCandidates[mapIndex].begin();
    multimap<int,ExtTrack>::iterator trkE = trackCandidates[mapIndex].end();
    //
    bool ADD = true;
    for( multimap<int,ExtTrack>::iterator trkI = trkB; trkI!=trkE; ++trkI){
      ExtTrack &trackC = (*trkI).second;
      bool SAME=true;
      for(int ip=0; ip<NEXTPLANES; ip++){
          bool isTrk = (ip < NEXTPLANES - _alg_nViewCal);
          if(          t_track.GetClusterX_ID(ip) != trackC.GetClusterX_ID(ip) )SAME=false;
          if(          t_track.GetClusterY_ID(ip) != trackC.GetClusterY_ID(ip) )SAME=false;
          if( isTrk && t_track.GetSensor(ip)      != trackC.GetSensor(ip)      )SAME=false;                   
      }
      if(SAME){
        ADD=false;
        break;
      }
    }
    //    if(_debug)cout <<" ADD ? "<<ADD<<endl;

      //      cout << "-----> INSERT  "<<endl;
      if(_debug){

          cout << endl;           
          cout<< " TRK chi2 "<<setw(13)<<t_track.chi2;
          cout <<" | al: ";
          for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
          cout <<" | ";
          cout << " X: ";
          //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.XGood(ip) ? "o" : "-");
          for(int ip=0; ip<t_track.nplanes; ip++){
              if(t_track.GetClusterX_ID(ip)>=0)
                  cout << setw(2)<<t_track.GetClusterX_ID(ip);
              else
                  cout << "__";
          }
          cout << " | ";
          cout << " Y: ";
          //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.YGood(ip) ? "o" : "-");
          for(int ip=0; ip<t_track.nplanes; ip++){
              if(t_track.GetClusterY_ID(ip)>=0)
                  cout << setw(2)<<t_track.GetClusterY_ID(ip);
              else
                  cout << "__";
          }
          cout << " S:";for(int ip=0; ip<t_track.nplanes; ip++)cout<< t_track.GetSensor(ip)+1;
          if(HBM)cout << " >> Merged";

      }


    if(ADD){

//      cout << " >> re-evaluate...";
      // ///////////////////////////////////////////////////////// FIT START
      int ifail=0;      
      t_track.ResetFit();  
      // t_track.Dump();
      t_track.Fit(0.,ifail,0);//
      if( ifail != 0 || TMath::IsNaN(t_track.chi2) ){
        // if(_debug){
        //   t_track.ResetFit();       
        //   t_track.Fit(0.,ifail,2);//
        // }
        continue;
      }
      
      // t_track.Fit(0.,ifail,0,1);//evaluate again taking into account the track angle
      // if( TMath::IsNaN(t_track.chi2) )continue;
      // if( ifail != 0 )continue;
      // ///////////////////////////////////////////////////////// FIT FINISH

      
      // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      // re-evaluated coordinates
      // and
      // evaluated other track info
      // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
      //      float VMS[22];VMS[0]=0.; // EM, 2245:13: warning: variable ?VMS? set but not used [-Wunused-but-set-variable]
      for(int ip=0; ip<NEXTPLANES; ip++){

        if(!t_track.XGood(ip)&&!t_track.YGood(ip))continue;
        
        int iclx = t_track.GetClusterX_ID(ip);
        int icly = t_track.GetClusterY_ID(ip);
        int ladder = t_track.GetLadder(ip);//ok
        int sensor = t_track.GetSensor(ip);// not yet assigned for calo

        bool isTrk = (ip < NEXTPLANES - _alg_nViewCal);   
        bool isXY = t_track.XGood(ip)*t_track.YGood(ip);
        bool isX  = t_track.XGood(ip) && !isXY;
        bool isY  = t_track.YGood(ip) && !isXY;
        
        //      cout <<endl<< ip << " xgood "<<t_track.xgood[ip]<<" ygood "<<t_track.ygood[ip];
        //      cout <<endl<< ip << " iclx "<<iclx<<" icly "<<icly;

        //      continue;
        //=========================================
        // evaluate angular factor
        //=========================================
        float factor = 0.;
        factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.axv[ip]/180.),2);
        factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.ayv[ip]/180.),2);
        factor += 1.;
        factor = TMath::Sqrt(factor);
        //      factor = 1.;///TEMP


        //=========================================
        // re-evaluate  coordinates
        //=========================================
        float xmABar[]={t_track.xv[ip],0.,0.,t_track.axv[ip],0.};//importante passare gli angoli e la coordinata!!
        float ymABar[]={t_track.yv[ip],0.,0.,t_track.ayv[ip],0.};//importante passare gli angoli e la coordinata!!
        float zmAB[]={t_track.zv[ip],0.,0.};
                
        if(isTrk){
            EvaluateClusterPosition_Tracker( iclx, icly, ladder, sensor, xmABar, ymABar, zmAB);           
        }else{      
            EvaluateClusterPosition_Calorimeter( (int)fmax((float)iclx,(float)icly), sensor, xmABar, ymABar, zmAB, t_track.al[4]);          
            // now the hiyt sensor is assigned
        }
        if(isXY){
          t_track.SetXY(ip,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
          t_track.SetZ(ip,zmAB[0]);
        }
        if(isX ){
          t_track.SetX( ip,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
          t_track.SetZ(ip,0.5*(zmAB[1]+zmAB[2]));
        }
        if(isY ){
          t_track.SetY( ip,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
          t_track.SetZ(ip,0.5*(zmAB[1]+zmAB[2]));
        }

        t_track.SetXGood(ip,(sensor == -1 ? 0 : iclx+1),ladder,sensor);
        t_track.SetYGood(ip,(sensor == -1 ? 0 : icly+1),ladder,sensor);
        
        //=========================================
        // re-evaluate  resolution (calorimeter only)
        // taking into account multiple scattering
        //=========================================
        // if(!isTrk && ip > NEXTPLANES - _alg_nViewCal){
        //   int nW = (int)((ip - 6 + 1)/2);
        //   float dW = 0.74*factor;//X0
        //   float dW_cm = 0.26*factor;//cm
        //   float dSi_cm = 0.; //cm
        //   if(isY) dSi_cm = ( ip%2 ? 0.228 : 0.428);
        //   dSi_cm *= factor;
        //   float theta = 0.;
        //   if(t_track.GetRigidity()!=0.)theta = 0.0136/t_track.GetRigidity();
        //   theta = theta * sqrt(dW) * (1.+0.038*log(dW));
        //   float Vt = theta*theta;            //cov(theta,theta)
        //   float Vy = Vt * dW_cm * dW_cm /3.; //cov(y,y)
        //   float Vty = 0.87 * Vt * Vy;        //cov(theta,y)
        //   float Vc = Vy + dSi_cm * dSi_cm * Vt + 2. * dSi_cm * Vty;
        //   float rms = sqrt( Vc + (nW>0 ? VMS[nW-1] : 0.) );    
        //   if(_debug)cout <<endl<< ip<<" nW "<<nW<<" res(MS) "<<rms;
        //   if(isY)t_track.resy[ip] = sqrt(rms*rms + t_track.resy[ip]*t_track.resy[ip]);
        //   else   t_track.resx[ip] = sqrt(rms*rms + t_track.resx[ip]*t_track.resx[ip]);
        //   if(isY)          VMS[nW] = rms*rms;
        // }

        //=========================================
        // evaluate other quantities
        //=========================================
        float smip,multmax;
        if(t_track.XGood(ip)){
          if(isTrk){
            TrkCluster *cluster = _trk_l1->GetCluster(iclx);
            float mip = TrkParams::GetMIP(cluster->GetLadder()-1,cluster->view-1);
            smip    = cluster->GetSignal()/(mip>0.?mip:1.); 
            multmax = cluster->GetMultiplicity()*10000+cluster->maxs;
          }else{
            ExtHit cluster = _cal_cl[iclx];
            smip    = cluster.signal;//mip
            multmax = cluster.mult*10000+lrint(cluster.coordPU);
          }
          t_track.dedx_x[ip]   = smip/(factor>0.?factor:1.); 
          t_track.multmaxx[ip] = multmax;
        }
        if(t_track.YGood(ip)){
          if(isTrk){
            TrkCluster *cluster = _trk_l1->GetCluster(icly);
            float mip = TrkParams::GetMIP(cluster->GetLadder()-1,cluster->view-1);
            smip    = cluster->GetSignal()/(mip>0.?mip:1.); 
            multmax = cluster->GetMultiplicity()*10000+cluster->maxs;
          }else{
            ExtHit cluster = _cal_cl[icly];
            smip    = cluster.signal;//mip
            multmax = cluster.mult*10000+lrint(cluster.coordPU);
          }
          t_track.dedx_y[ip]   = smip/(factor>0.?factor:1.); 
          t_track.multmaxy[ip] = multmax;
        }
      }
      ///////////////////////////////////////////////////////// end additional info
 

      // ///////////////////////////////////////////////////////// RE-FIT START
      t_track.Fit(0.,ifail,0);//evaluate again taking into account the track angle
      if( TMath::IsNaN(t_track.chi2) )continue;
      if( ifail != 0 )continue;
      // ///////////////////////////////////////////////////////// RE-FIT FINISH
      if(_debug){
//      cout << endl;
//              cout << " >> al: ";for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
//              cout << " chi2: "<<setw(10)<< t_track.chi2;
        //      cout << endl;
          cout << endl;
//                cout << " >> al: ";for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
//                cout << " chi2: "<<setw(10)<< t_track.chi2;
//                //            cout << endl;
          cout<< " TRK chi2 "<<setw(13)<<t_track.chi2;
          cout <<" | al: ";
          for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
          cout <<" | ";
          cout << " X: ";
          //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.XGood(ip) ? "o" : "-");
          for(int ip=0; ip<t_track.nplanes; ip++){
              if(t_track.GetClusterX_ID(ip)>=0)
                  cout << setw(2)<<t_track.GetClusterX_ID(ip);
              else
                  cout << "__";
          }
          cout << " | ";
          cout << " Y: ";
          //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.YGood(ip) ? "o" : "-");
          for(int ip=0; ip<t_track.nplanes; ip++){
              if(t_track.GetClusterY_ID(ip)>=0)
                  cout << setw(2)<<t_track.GetClusterY_ID(ip);
              else
                  cout << "__";
          }
          cout << " S:";for(int ip=0; ip<t_track.nplanes; ip++)cout<< t_track.GetSensor(ip)+1;
          cout << " >> ADD ";
      }

     _caloTj->DoTrack(t_track.al,t_track.zini);
      for(int ip=0; ip<_caloTj->npoint; ip++){
        t_track.xGF[ip] = _caloTj->x[ip];
        t_track.yGF[ip] = _caloTj->y[ip];
      }

      //              cout << "ADD "<<endl;
      trackCandidates[mapIndex].insert(pair<int,ExtTrack>(t_track.GetNX()+t_track.GetNY(), ExtTrack(t_track)));
      if(t_track.GetNX()+t_track.GetNY()>5 && t_track.chi2<chi2Min) chi2Min=t_track.chi2;
      iCand++;
    }//end add-candidate condition
            
  }//end first loop on candidates
  if(_debug)cout <<endl<< "n.candidates "<<trackCandidates[mapIndex].size()<<endl;

//  cout << "CHI2 min "<<chi2Min<<endl;
  //=================================================================
  // ----------------------------------------------------------------
  // chi2 selection
  // ----------------------------------------------------------------
  //=================================================================
  pair <multimap<int,ExtTrack>::iterator, multimap<int,ExtTrack>::iterator> nxny[NEXTVIEWS];
//  for(int inxny=0; inxny<NEXTVIEWS; inxny++){
  bool NMTI =false;//not-minimal track inserted
  for(int inxny=NEXTVIEWS-1; inxny>=nClFixX+nClFixY; inxny--){
    nxny[inxny] = trackCandidates[mapIndex].equal_range(inxny);
//    cout << " --- n "<<inxny<<endl;
    for ( multimap<int,ExtTrack>::iterator it = nxny[inxny].first; it!=nxny[inxny].second; ++it){
      ExtTrack &trackC = (*it).second;
//      cout << " NX+NY "<< trackC.GetNX()+trackC.GetNY() <<" chi2 "<<trackC.chi2<< " R(GV) "<<trackC.GetRigidity()<<endl;
      //            if( trackC.chi2==0 )continue;
      if( (trackC.GetNX()+trackC.GetNY() > 5 && trackC.chi2>chi2Min  ) ||
          (trackC.GetNX()+trackC.GetNY() == 5 && NMTI ) ||
          false)continue;      

      trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(inxny, ExtTrack(trackC)));
//      cout << " insert "<<endl;

      if( trackC.GetNX()+trackC.GetNY() > 5 )NMTI=true;
    }
  }
  mapIndex++;

  //  if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;

  if(trackCandidates[mapIndex].size()==0)return;//no good candidates

  if(_debug){     
    
    cout<< " Selected-track SUMMARY:"<<endl;

    trkBegin = trackCandidates[mapIndex].begin();
    trkEnd   = trackCandidates[mapIndex].end();
    for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
      ExtTrack &trackCand = (*trkIter).second;
      cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
      cout <<" | al: ";
      for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
      cout <<" | ";
      cout << " X: ";
      //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
      for(int ip=0; ip<trackCand.nplanes; ip++){
        if(trackCand.GetClusterX_ID(ip)>=0)
          cout << setw(2)<<trackCand.GetClusterX_ID(ip);
        else
          cout << "__";
      }
      cout << " | ";
      cout << " Y: ";
      //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
      for(int ip=0; ip<trackCand.nplanes; ip++){
        if(trackCand.GetClusterY_ID(ip)>=0)
          cout << setw(2)<<trackCand.GetClusterY_ID(ip);
        else
          cout << "__";
      }
      cout << " S:";for(int ip=0; ip<trackCand.nplanes; ip++)cout<< trackCand.GetSensor(ip)+1;
      cout <<"  MDR="<<1./sqrt(trackCand.coval[4][4])<<endl;
      cout << endl;
      // for(int ip=0; ip<trackCand.nplanes; ip++)
      //        if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";       cout << endl;

    }
  }
  //=================================================================
  // ----------------------------------------------------------------
  // save tracks
  // ----------------------------------------------------------------
  //=================================================================
  trkBegin = trackCandidates[mapIndex].begin();
  trkEnd   = trackCandidates[mapIndex].end();
  iCand = 0;
  // TClonesArray &tracks = *_trkArray;
  // tracks.Clear();
  for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){            
    ExtTrack trackCand = ((*trkIter).second);
    //    trackCand.seqno = iCand;
    //    trackCand.Dump();
    new(tracks[iCand]) ExtTrack(trackCand);
    iCand++;
  };


};
/**
 * 
 */

void ExtTrkingAlg::ProcessEvent2(Bool_t force){


//    _debug=true;

  if(_debug)cout << " |---------------------------------------------------| "<<endl;
  if(_debug)cout << " void ExtTrkingAlg::ProcessEvent2("<<force<<") "<<endl;

  
  if(_debug && !_trk_l1)cout << " Missing TrkLevel1 object --- ciao! "<<endl;
  if(!_trk_l1)return;


  //===========================
  // pre-selection  
  //===========================
  
  Bool_t RETRACK = CheckEvent();

  if(_debug)cout <<" Pre-selection result: "<<RETRACK<<endl;
  if( !force && !RETRACK )return;

  NEXTVIEWS  = 12 +_alg_nViewCal;
  NEXTPLANES = 6  +_alg_nViewCal;

  //===========================
  // tracking
  //===========================

  // TClonesArray &tracks = *_trkArray;
  // tracks.Clear();
  TClonesArray &tracks = *_trkArray;
  tracks.Clear("C");

  int maxSize = (5+5)*(_alg_nViewCal+1);//boh?!
  multimap<int,ExtTrack> trackCandidates[maxSize]; // map: last hit view vs TrkTrack
  int mapIndex = 0;
  multimap<int,ExtTrack>::iterator trkBegin;
  multimap<int,ExtTrack>::iterator trkEnd;
  int iCand;

  // ------------------------------------------------
  // fill a map VIEW-CLUSTER
  // ------------------------------------------------
  multimap<int,int> clusterMap;
  FillClusterMap(clusterMap,_trk_l1,0);// trk clusters    
  FillClusterMap(clusterMap,_cal_l1,12);// calorimeter clusters
  

  // ------------------------------------------------
  // define iterators over clusters on the same view
  // ------------------------------------------------
  pair <multimap<int,int>::iterator, multimap<int,int>::iterator> ret[NEXTVIEWS];
  multimap<int,int>::iterator it[NEXTVIEWS];
  for(int iv=0; iv<NEXTVIEWS; iv++)ret[iv] = clusterMap.equal_range(iv);


  //====================================================================
  // PARAMETERS
  //====================================================================
  // ------------------------------------------------
  // number of points per candidate
  // ------------------------------------------------
  int nClFixX =  _alg_nClFixX; // n.hits required on X view
  int nClFixY = _alg_nClFixY; // n.hits required on X view
  int nTrackMAX =  _alg_nTrackMAX; // 


  /////////////////////////////////////////////////////////////////////////
  //
  //
  // start loop over n.calo planes to be included
  //
  /////////////////////////////////////////////////////////////////////////

  for(int nViewCal=0; nViewCal<=_alg_nViewCal; nViewCal++){

      NEXTVIEWS  = 12 +nViewCal;
      NEXTPLANES = 6  +nViewCal;
      
      if(_debug)cout<<">>>>>>>>>>> INCLUDE "<<nViewCal<<" CALORIMETER VIEWS "<<endl;
//      cout<<">>>>>>>>>>> INCLUDE "<<nViewCal<<" CALORIMETER VIEWS "<<endl;
      //==================================================
      // ------------------------------------------------
      // start: one candidate for each cluster
      // ------------------------------------------------
      //==================================================
      if(_debug)cout<<"-------- MINIMAL TRACK-CANDIDATES --------";
      if(_debug)cout<<nClFixX<<"X"<<nClFixY<<"Y"<<endl;
      

      for(int iv=0; iv<NEXTVIEWS; iv++){                                      //loop over all views
          
          if( !(iv%2) && (int)(0.5*iv) >  NEXTVIEWS/2-nClFixX)continue;
          if(  (iv%2) && (int)(0.5*iv) >  NEXTVIEWS/2-nClFixY)continue;
          
          if(iv >= NEXTVIEWS - nViewCal )continue; //skip calo hits 
          
          for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view
              
              
              int icl = it[iv]->second; 

              TrkCluster *cluster = _trk_l1->GetCluster(icl);
              //      cluster->Dump();
              int ladder = cluster->GetLadder()-1; //
              int plane = cluster->GetPlane()-1;
              bool isY = (cluster->view)%2;// F77 1-12 ??
              bool isX = !isY;
      

              /// TRACKER ==============================================================================
              for (int is=0; is<2; is++){                                  //loop over sensors

                  //cout <<endl<< " view "<<iv<<" plane "<<plane<<" cluster "<<icl<<" on sensor "<<is<<endl;
        
                  ExtTrack& t_track =  *_extTrack;

                  t_track.ResetXY();
//                for(int ip =0 ; ip<t_track.nplanes; ip++ )t_track.SetZ(ip,_zMech[ip]);
//                t_track.Dump();


                  // >>>> insert point
                  float xmABar[]={0.,0.,0.,0.,0.}; // inizializzare a zero l'angolo, come prima stima per il PFA
                  float ymABar[]={0.,0.,0.,0.,0.};
                  float zmAB[]={t_track.zm[plane],0.,0.};//mechanical z-coordinate
                  EvaluateClusterPosition_Tracker( (isX?icl:-1), (isY?icl:-1), ladder, is, xmABar, ymABar, zmAB); 
        
                  if( isX ){
                      t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
                      t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
                      t_track.SetXGood(plane,icl+1,ladder,is);  
                      //          t_track.SetYGood(plane,0,ladder,is);  
                  }
                  if( isY ){
                      t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
                      t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
                      t_track.SetYGood(plane,icl+1,ladder,is);
                      //          t_track.SetXGood(plane,0,ladder,is);
                  }
                  trackCandidates[mapIndex].insert(pair<int,ExtTrack>(iv, t_track));

              }
      

          }
      }


      //==============================================================
      // -------------------------------------------------------------
      // next: add views to each candidate, up to nClFixX+nClFixY
      // -------------------------------------------------------------
      //==============================================================
//      for( ;mapIndex < nClFixX+nClFixY-1 ;++mapIndex){
      for( int iii=0; iii< nClFixX+nClFixY-1; iii++){

//        if(mapIndex>0)trackCandidates[mapIndex-1].clear();

          trkBegin = trackCandidates[mapIndex].begin();
          trkEnd   = trackCandidates[mapIndex].end();
          if(_debug)cout<<"MAP "<<mapIndex<<" ----------- n.candidates "<<trackCandidates[mapIndex].size()<<endl;

          for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){


              ExtTrack trackCand = (*trkIter).second;
              int lastView = (*trkIter).first; //last inserted view
              bool lastIsTrk = (lastView < NEXTVIEWS - nViewCal);
              int lastPlane = (int)(0.5*lastView);
              if(!lastIsTrk)lastPlane =lastView-6;
              int lastLadder = trackCand.GetLadder(lastPlane);
              int lastSensor = trackCand.GetSensor(lastPlane);


//            cout <<" >>> lastView "<<lastView<< " NEXTVIEWS "<<NEXTVIEWS<<endl;

              for(int iv=lastView+1; iv<NEXTVIEWS; iv++){                             //loop over next views

                  bool isTrk = (iv < NEXTVIEWS - nViewCal);
                  bool isX = (iv)%2;
                  bool isY = !isX;
        
                     
                  // if(  (iv%2) && (int)(0.5*iv) > 6-nClFixX+trackCand.GetNX())continue; //not enough x-view left
                  // if( !(iv%2) && (int)(0.5*iv) > 6-nClFixY+trackCand.GetNY())continue; //not enough y-view left

                  // if(  (iv%2) && trackCand.GetNX()==nClFixX )continue;//ok?
                  // if( !(iv%2) && trackCand.GetNY()==nClFixY )continue;//ok?

                  if( isX && trackCand.GetNX()==nClFixX )continue;//ok?
                  if( isY && trackCand.GetNY()==nClFixY )continue;//ok?


                  if( isX && (int)(0.5*iv) >  6+nViewCal/2     - nClFixX+trackCand.GetNX())continue; //not enough x-view left
                  if( isY && (int)(0.5*iv) >  6+(nViewCal+1)/2 - nClFixY+trackCand.GetNY())continue; //not enough y-view left

                  for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view

                      int icl = it[iv]->second;

//                    cout <<" >>> icl "<<icl<<endl;

                      if(isTrk){
                          /// TRACKER ==================================================
            

                          TrkCluster *cluster = _trk_l1->GetCluster(icl);
                          int ladder = cluster->GetLadder()-1;//(_trk_l1->GetCluster(icl)->maxs-1)%1024 + 1;
                          int plane = cluster->GetPlane()-1;
                          // bool isY = (cluster->view)%2;
                          // bool isX = !isY;
            
                          if( plane==lastPlane && ladder!=lastLadder)continue;
            

                          for (int is=0; is<2; is++){                                  //loop over sensors
                  
                              if( plane==lastPlane && is!=lastSensor)continue;
              
                              // if(_debug)cout<<"MAP "<<mapIndex<<" last v "<<lastView;
                              // if(_debug)cout<<" +cl "<<icl<<" s"<<is<<endl;
              
                              ExtTrack t_track = ExtTrack();
                              trackCand.Copy(t_track);                             //copy the candidate parameters...
              
                              bool isXY = (isX && t_track.YGood(plane)) || (isY && t_track.XGood(plane)) ;
              

                              float xmABar[]={0.,0.,0.,0.,0.}; //importante inizializzare a zero, come prima stima per il PFA
                              float ymABar[]={0.,0.,0.,0.,0.};
                              float zmAB[]={t_track.zm[plane],0.,0.};

                              if( isX ){
                                  int iclx = icl;
                                  int icly = (isXY ? t_track.GetClusterY_ID(plane) : -1 );

                                  EvaluateClusterPosition_Tracker( iclx, icly , ladder, is, xmABar, ymABar, zmAB); 

                                  if(isXY){
                                      t_track.SetXY(plane,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
                                      t_track.SetZ(plane,zmAB[0]);
                                  }else{
                                      t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
                                      t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
                                  }
                                  t_track.SetXGood(plane,iclx+1,ladder,is);         //...add a point...
                                  t_track.SetYGood(plane,icly+1,ladder,is);
                              }
                              //
                              if( isY ){
                                  int iclx = (isXY ? t_track.GetClusterX_ID(plane) : -1 );
                                  int icly = icl;

                                  EvaluateClusterPosition_Tracker( iclx, icly , ladder, is, xmABar, ymABar, zmAB); 

                                  if(isXY){
                                      t_track.SetXY(plane,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
                                      t_track.SetZ(plane,zmAB[0]);
                                  }else{
                                      t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
                                      t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
                                  };
                                  t_track.SetYGood(plane,icly+1,ladder,is);
                                  t_track.SetXGood(plane,iclx+1,ladder,is);
                              }

                              trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(iv, t_track)); //...and store a new candidate
                          };//end loop over sensors
            
                      }else{          
                          /// CALORIMETER =================================================
                          ExtHit cluster = _cal_cl[icl];
                          int plane = 6 + cluster.view;

                          if( plane==lastPlane )continue;

                          int ladder = (int)(cluster.coordPU / 32);
                          int sensor = -1; //not yet known
                          bool isX = (cluster.view)%2;
                          bool isY = !isX;

                          ExtTrack t_track = ExtTrack();
                          trackCand.Copy(t_track);   

                          float xmABar[]={0.,0.,0.,0.,0.}; //importante inizializzare a zero, come prima stima per il PFA
                          float ymABar[]={0.,0.,0.,0.,0.};
                          float zmAB[]={0.,0.,0.};
                          EvaluateClusterPosition_Calorimeter( icl, sensor, xmABar, ymABar, zmAB); 


                          if( isX ){
                              t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
                              t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
                              t_track.SetXGood(plane,icl+1, ladder, sensor);         //...add a point.
                              t_track.SetYGood(plane,    0, ladder, sensor);     //...add a point...
                          }        
                          if( isY ){
                              t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
                              t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
                              t_track.SetYGood(plane,icl+1, ladder, sensor);
                              t_track.SetXGood(plane,    0, ladder, sensor);
                          }
                          trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(iv, t_track)); //...and store a new ca

                      }

                      if( trackCandidates[mapIndex+1].size() > (uint)nTrackMAX ){ //EM, 2826:50: warning: comparison between signed and unsigned integer expressions [-Wsign-compare]
                          if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex+1].size()<<" > "<<nTrackMAX<<endl;
                          return;//to many candidates
                      }

                  };//end loop over clusters on the view iv
              };//endl loop over views
          }//end loop over candidates

          mapIndex++;

      }//end iterations
  
      if(_debug)cout<<"end of cluster inclusion"<<endl;
      if(_debug)cout<<"MAP "<<mapIndex<<" ----------- n.candidates "<<trackCandidates[mapIndex].size()<<endl;
//      if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;


      //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
      if(trackCandidates[mapIndex].size()==0){
//        mapIndex++;
          continue;//increment one calo plane
      }
      //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+


      //===============================================================================================
      // -------------------------------------------------------------
      // first selection of track candidates, based on chi2
      // -------------------------------------------------------------
      //===============================================================================================
      if(_debug)cout<<"-------- CHI2 --------"<<endl;
//      if(mapIndex>0)trackCandidates[mapIndex-1].clear();
      trkBegin = trackCandidates[mapIndex].begin();
      trkEnd   = trackCandidates[mapIndex].end();
      mapIndex++;
      iCand = 0;
      for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
          int lastView = (*trkIter).first;
          ExtTrack &trackCand = (*trkIter).second;        
          int ifail=0;
          

          trackCand.ResetFit();
          trackCand.Fit(0.,ifail,0);
          if(ifail!=0)trackCand.ResetFit();
          
//        trackCand.Dump();


          // -----------------------------------
          // first selection of track candidates
          // -----------------------------------
          if(TMath::IsNaN(trackCand.chi2))continue;
          if(trackCand.chi2 <= 0. || trackCand.chi2 > 1.e7)continue; 
          //    if(trackCand.chi2 < 0. || trackCand.chi2 > 10. )continue; 
          iCand++;    
          
          if(_debug){     
              cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
              cout <<" | al: ";
              for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
              cout <<" | ";
              cout << " X: ";
              for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
              // for(int ip=0; ip<trackCand.nplanes; ip++)
              // if(trackCand.GetClusterX_ID(ip)>=0)cout << trackCand.GetClusterX_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";
              cout << " | ";
              cout << " Y: ";
              for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
              cout << endl;
              // for(int ip=0; ip<trackCand.nplanes; ip++)
              //        if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";       cout << endl;
          }
          trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, trackCand));

      }


      //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
      if(trackCandidates[mapIndex].size()==0){
          mapIndex++;
          continue;//increment one calo plane
      }
      //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

      if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;


      //=================================================================
      // ----------------------------------------------------------------
      // CALORIMETER
      // ----------------------------------------------------------------
      //=================================================================
      if(_cal_l1){
          if(_debug)cout<<"-------- CALORIMETER --------"<<endl;

          CaloLevel1* l1=_cal_l1;//event->GetCaloLevel1();
          //    if(!l1)return;
          vector<float> caloChi2;
          vector<int> caloHit; int caloHitMAX=0;
          vector<float> caloMip;
//        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
          trkBegin = trackCandidates[mapIndex].begin();
          trkEnd   = trackCandidates[mapIndex].end();
          for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){          
              ExtTrack trackCand = ((*trkIter).second);   
            
              // trackCand.DoTrack(&caloTj);//integrate the trajectory through the calorimeter
              _caloTj->DoTrack(trackCand.al,trackCand.zini);
              //      _caloTj->Dump();
      
              float chi2 = 0.;
              int  nhit = 0;
              float cmip = 0.;
              for(Int_t ih=0; ih<l1->istrip; ih++){
                  Int_t view = -1;
                  Int_t plane = -1;
                  Int_t strip = -1;
                  Float_t mip = l1->DecodeEstrip(ih,view,plane,strip);
                  if(view<0)continue;
                  if(plane<0)continue;
                  if(strip<0)continue;
                  if(mip<=0)continue;
                  float dxy = _caloCoord[view][plane][strip]-(view?_caloTj->y[plane*2+(view?0:1)]:_caloTj->x[plane*2+(view?0:1)]);
                  if(plane<11) chi2+= mip * TMath::Power(dxy,2.);
                  //          caloZ[plane*2+(view?0:1)]=st.GetZ();
                  if( fabs(dxy) < 2. && plane<11)nhit++;
                  if( fabs(dxy) < 2. && plane<11)cmip+=mip;
              };
              if(l1->istrip>0)chi2 = chi2/l1->istrip;
            
              caloHit.push_back(nhit); 
              if(nhit>caloHitMAX)caloHitMAX=nhit;

              if(_debug){           
                  cout << " CALO ";
                  //    cout << " chi2 "<< chi2 <<" ";
                  cout << " nhit (< 11th plane) "<< nhit <<" ";
                  //    cout << " mip  "<< cmip <<" ";
                  cout <<endl;
              };
          }

          //=================================================================
          // ----------------------------------------------------------------
          // selection of track candidates, based on n.calorimeter hits
          // ----------------------------------------------------------------
          //=================================================================
//        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
          trkBegin = trackCandidates[mapIndex].begin();
          trkEnd   = trackCandidates[mapIndex].end();
          mapIndex++;
          iCand = 0;
          for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
              int lastView = (*trkIter).first;
              ExtTrack &trackCand = (*trkIter).second;
              if(caloHit[iCand]==caloHitMAX)trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, trackCand));
              iCand++;
          }

          //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
          if(trackCandidates[mapIndex].size()==0){
              mapIndex++;
              continue;//increment one calo plane
          }
          //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

          if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;

      }// end if calo event


      //=================================================================
      // ----------------------------------------------------------------
      // TOF
      // ----------------------------------------------------------------
      //=================================================================
      if(_tof_l2){
          if(_debug)cout<<"-------- TOF --------"<<endl;
//        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
          trkBegin = trackCandidates[mapIndex].begin();
          trkEnd   = trackCandidates[mapIndex].end();
          ToFLevel2* tl2=_tof_l2;//event->GetToFLevel2(); 
          //   if(!tl2)return;
          iCand =0;
          vector<int> tofHit; int tofHitMAX=0; 
          for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
              //            int lastView = (*trkIter).first;
              ExtTrack &trackCand = (*trkIter).second;

              _tgf->DoTrack(trackCand.al);//<<<< integrate the trajectory

              for(int ip=0; ip<TrkParams::nGF; ip++){
                  trackCand.xGF[ip] = _tgf->x[ip];
                  trackCand.yGF[ip] = _tgf->y[ip];
              }

              int nhit =0;
              for(int iToF=0; iToF<6; iToF++){
              
                  int iGF = iToF;
                  if(iToF>3)iGF = iToF+8;
                  int iPaddle = tl2->GetPaddleIdOfTrack( trackCand.xGF[iGF], trackCand.yGF[iGF], iToF);
                  if( tl2->HitPaddle(iToF,iPaddle) )nhit++;
              
              }
              iCand++;
            
              tofHit.push_back(nhit); 
              if(nhit>tofHitMAX)tofHitMAX=nhit;

              if(_debug){
                  cout << " TOF nhit "<< nhit <<" ";
                  cout <<endl;
              }
            
    
          }
          //=================================================================
          // ----------------------------------------------------------------
          // selection of track candidates, based on n.tof hits
          // ----------------------------------------------------------------
          //=================================================================
//        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
          trkBegin = trackCandidates[mapIndex].begin();
          trkEnd   = trackCandidates[mapIndex].end();
          mapIndex++;
          iCand = 0;
          for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
              int lastView = (*trkIter).first;
              ExtTrack &trackCand = (*trkIter).second;
              if(tofHit[iCand]==tofHitMAX)trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, ExtTrack(trackCand)));
              iCand++;

          }
          if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;

          //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
          if(trackCandidates[mapIndex].size()==0){
              mapIndex++;
              continue;//increment one calo plane
          }
          //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

      }

      if(_debug){     
    
          cout<< " Minimal-track SUMMARY:"<<endl;

//        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
          trkBegin = trackCandidates[mapIndex].begin();
          trkEnd   = trackCandidates[mapIndex].end();
          for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
              ExtTrack &trackCand = (*trkIter).second;
              cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
              cout <<" | al: ";
              for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
              cout <<" | ";
              cout << " X: ";
              //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
              for(int ip=0; ip<trackCand.nplanes; ip++){
                  if(trackCand.GetClusterX_ID(ip)>=0)
                      cout << setw(2)<<trackCand.GetClusterX_ID(ip);
                  else
                      cout << "__";
              }
              cout << " | ";
              cout << " Y: ";
              //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
              for(int ip=0; ip<trackCand.nplanes; ip++){
                  if(trackCand.GetClusterY_ID(ip)>=0)
                      cout << setw(2)<<trackCand.GetClusterY_ID(ip);
                  else
                      cout << "__";
              }
              cout << " S:";for(int ip=0; ip<trackCand.nplanes; ip++)cout<< trackCand.GetSensor(ip)+1;
              cout << endl;

          }
      }
      //=================================================================
      // ----------------------------------------------------------------
      // TRACK MERGING
      // ----------------------------------------------------------------
      //=================================================================
      if(_debug)cout<<"-------- MERGING --------";
//        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
      trkBegin = trackCandidates[mapIndex].begin();
      trkEnd   = trackCandidates[mapIndex].end();
      mapIndex++;
      float chi2Min=1e10;
      iCand = 0;
      for( multimap<int,ExtTrack>::iterator trkIter1 = trkBegin; trkIter1!=trkEnd; ++trkIter1){ //loop over tracks (1)
    
          ExtTrack &trackCand1 = (*trkIter1).second;// get candidate 1...
            
          ExtTrack t_track = ExtTrack();// ...create a new ExtTrack instance...
          trackCand1.Copy(t_track);     // ...and make a copy
                   
          bool HBM = false; //has been merged

          for( multimap<int,ExtTrack>::iterator trkIter2 = trkBegin; trkIter2!=trkEnd; ++trkIter2){//loop over tracks (2)
              
              ExtTrack &trackCand2 = (*trkIter2).second; // get candidate 2...

              bool CBM = true;
              for(int ip=0; ip<NEXTPLANES; ip++){ //loop over planes
                  if( t_track.XGood(ip)  &&  
                      trackCand2.XGood(ip) &&  
                      ( t_track.GetClusterX_ID(ip)!=trackCand2.GetClusterX_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
                      true )CBM=false; // if both have a x-hit and it is not the same ---> they cannot be merged
                  if(!CBM)break;       
                  if( t_track.YGood(ip)  && 
                      trackCand2.YGood(ip) && 
                      ( t_track.GetClusterY_ID(ip)!=trackCand2.GetClusterY_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
                      true )CBM=false; // if both have a y-hit and it is not the same ---> they cannot be merged
                  if(!CBM)break;        
              }
              //      if(_debug)cout << " can be merged ? "<<CBM<<endl;
              if(!CBM)continue; //go to the next candidate
              ///////////
              //merging//
              ///////////    
              // cout <<endl<<" check merging al (1) ";
              // for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
              // cout <<endl<<"       with    al (2) ";
              // for(int i=0; i<5; i++)cout <<setw(10)<< trackCand2.al[i];
              for(int ip=0; ip<NEXTPLANES; ip++){             
                  if( !(t_track.XGood(ip)) &&  trackCand2.XGood(ip) ){            
                      int ygood_temp = t_track.ygood[ip];
          
                      if(t_track.YGood(ip)){ 
                          t_track.SetXY(ip,trackCand2.xm[ip],t_track.ym[ip],trackCand2.resx[ip],t_track.resy[ip]);
                          //        cout<<"confronta "<<trackCand2.zm[ip]<<" - "<<t_track.zm[ip]<<" = "<<trackCand2.zm[ip]-t_track.zm[ip]<<endl;
                      }else{ 
                          t_track.SetX(ip,trackCand2.xma[ip],trackCand2.xmb[ip],trackCand2.yma[ip],trackCand2.ymb[ip],trackCand2.resx[ip]);
                          t_track.SetZ(ip,0.5*(trackCand2.zma[ip]+trackCand2.zmb[ip]));
                      }
                      t_track.xgood[ip] = trackCand2.xgood[ip];//assign cluster + sensor
                      t_track.ygood[ip] = ygood_temp;//assign cluster + sensor

                      HBM = true;        
                  }             
                  if( !(t_track.YGood(ip)) &&  trackCand2.YGood(ip) ) {           

                      int xgood_temp = t_track.xgood[ip];
          
                      if(t_track.XGood(ip)){ 
                          t_track.SetXY(ip,t_track.xm[ip],trackCand2.ym[ip],t_track.resx[ip],trackCand2.resy[ip]);       
                          //        cout<<"confronta "<<trackCand2.zm[ip]<<" - "<<t_track.zm[ip]<<" = "<<trackCand2.zm[ip]-t_track.zm[ip]<<endl;
                      }else{ 
                          t_track.SetY(ip,trackCand2.xma[ip],trackCand2.xmb[ip],trackCand2.yma[ip],trackCand2.ymb[ip],trackCand2.resy[ip]);
                          t_track.SetZ(ip,0.5*(trackCand2.zma[ip]+trackCand2.zmb[ip]));
                      }
                      t_track.ygood[ip] = trackCand2.ygood[ip];//assign cluster + sensor 
                      t_track.xgood[ip] = xgood_temp;//assign cluster + sensor
                      HBM = true;
                  }
              };         
      
          }//end loop over tracks (2)
          //
          //check if track has been already inserted
          //
//        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
          multimap<int,ExtTrack>::iterator trkB = trackCandidates[mapIndex].begin();
          multimap<int,ExtTrack>::iterator trkE = trackCandidates[mapIndex].end();
          //
          bool ADD = true;
          for( multimap<int,ExtTrack>::iterator trkI = trkB; trkI!=trkE; ++trkI){
              ExtTrack &trackC = (*trkI).second;
              bool SAME=true;
              for(int ip=0; ip<NEXTPLANES; ip++){
                  bool isTrk = (ip < NEXTPLANES - nViewCal);
                  if(          t_track.GetClusterX_ID(ip) != trackC.GetClusterX_ID(ip) )SAME=false;
                  if(          t_track.GetClusterY_ID(ip) != trackC.GetClusterY_ID(ip) )SAME=false;
                  if( isTrk && t_track.GetSensor(ip)      != trackC.GetSensor(ip)      )SAME=false;                   
              }
              if(SAME){
                  ADD=false;
                  break;
              }
          }
          //    if(_debug)cout <<" ADD ? "<<ADD<<endl;

              //      cout << "-----> INSERT  "<<endl;
          if(_debug){

              cout << endl;               
//            cout<< " TRK chi2 "<<setw(13)<<t_track.chi2;
              cout<< " TRK chi2 "<<".............";
              cout <<" | al: ";
//            for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
              for(int i=0; i<5; i++)cout <<"..........";
              cout <<" | ";
              cout << " X: ";
              //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.XGood(ip) ? "o" : "-");
              for(int ip=0; ip<t_track.nplanes; ip++){
                  if(t_track.GetClusterX_ID(ip)>=0)
                      cout << setw(2)<<t_track.GetClusterX_ID(ip);
                  else
                      cout << "__";
              }
              cout << " | ";
              cout << " Y: ";
              //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.YGood(ip) ? "o" : "-");
              for(int ip=0; ip<t_track.nplanes; ip++){
                  if(t_track.GetClusterY_ID(ip)>=0)
                      cout << setw(2)<<t_track.GetClusterY_ID(ip);
                  else
                      cout << "__";
              }
              cout << " S:";for(int ip=0; ip<t_track.nplanes; ip++)cout<< t_track.GetSensor(ip)+1;
              if(HBM)cout << " >> Merged";
//            cout << endl;
          }       

          if(ADD){
//            cout << " >> re-evaluate...";
//                cout << " >> Merged? "<<HBM;
//                cout << "| X:";
//                for(int ip=0; ip<NEXTPLANES; ip++){
//                    if(t_track.GetClusterX_ID(ip)>=0)
//                        cout << setw(2)<<t_track.GetClusterX_ID(ip);
//                    else
//                        cout << "__";
//                }
//                //    for(int ip=0; ip<NEXTPLANES; ip++)cout<< t_track.XGood(ip);
//                cout << "| Y:";
//                for(int ip=0; ip<NEXTPLANES; ip++){
//                    if(t_track.GetClusterY_ID(ip)>=0)
//                        cout << setw(2)<<t_track.GetClusterY_ID(ip);
//                    else
//                        cout << "__";
//                }
//                //    for(int ip=0; ip<NEXTPLANES; ip++)cout<< t_track.YGood(ip);     
//                cout << " S:";for(int ip=0; ip<NEXTPLANES; ip++)cout<< t_track.GetSensor(ip)+1;
//                cout << endl;
              
              // ///////////////////////////////////////////////////////// FIT START
              int ifail=0;      
              t_track.ResetFit();  
              // t_track.Dump();
              t_track.Fit(0.,ifail,0);//
              if( ifail != 0 || TMath::IsNaN(t_track.chi2) ){
                  // if(_debug){
                  //   t_track.ResetFit();       
                  //   t_track.Fit(0.,ifail,2);//
                  // }
                  continue;
              }
      
              // t_track.Fit(0.,ifail,0,1);//evaluate again taking into account the track angle
              // if( TMath::IsNaN(t_track.chi2) )continue;
              // if( ifail != 0 )continue;
              // ///////////////////////////////////////////////////////// FIT FINISH

      
              // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
              // re-evaluated coordinates
              // and
              // evaluated other track info
              // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
              //float VMS[22];VMS[0]=0.;  //EM, 3295:14: warning: variable ?VMS? set but not used [-Wunused-but-set-variable]
              for(int ip=0; ip<NEXTPLANES; ip++){

                  if(!t_track.XGood(ip)&&!t_track.YGood(ip))continue;
        
                  int iclx = t_track.GetClusterX_ID(ip);
                  int icly = t_track.GetClusterY_ID(ip);
                  int ladder = t_track.GetLadder(ip);//ok
                  int sensor = t_track.GetSensor(ip);//not assigned yet for calo

                  bool isTrk = (ip < NEXTPLANES - nViewCal);      
                  bool isXY = t_track.XGood(ip)*t_track.YGood(ip);
                  bool isX  = t_track.XGood(ip) && !isXY;
                  bool isY  = t_track.YGood(ip) && !isXY;
        
                  //    cout <<endl<< ip << " xgood "<<t_track.xgood[ip]<<" ygood "<<t_track.ygood[ip];
                  //    cout <<endl<< ip << " iclx "<<iclx<<" icly "<<icly;

                  //    continue;
                  //=========================================
                  // evaluate angular factor
                  //=========================================
                  float factor = 0.;
                  factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.axv[ip]/180.),2);
                  factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.ayv[ip]/180.),2);
                  factor += 1.;
                  factor = TMath::Sqrt(factor);
                  //    factor = 1.;///TEMP

                  //=========================================
                  // re-evaluate  coordinates
                  //=========================================
                  float xmABar[]={t_track.xv[ip],0.,0.,t_track.axv[ip],0.};//importante passare gli angoli!!
                  float ymABar[]={t_track.yv[ip],0.,0.,t_track.ayv[ip],0.};//importante passare gli angoli!!
                  float zmAB[]={t_track.zv[ip],0.,0.};
                
                  if(isTrk){
                      EvaluateClusterPosition_Tracker( iclx, icly, ladder, sensor, xmABar, ymABar, zmAB);         
                  }else{
//                    cout <<" ------------ EvaluateClusterPosition_Calorimeter ----------------"<<endl;
//                    cout << "IN: "<<(int)fmax((float)iclx,(float)icly)<<endl;
//                    cout << "IN: "<<sensor<<endl;
//                    cout << "IN: "<<xmABar[0]<<endl;
//                    cout << "IN: "<<ymABar[0]<<endl;
//                    cout << "IN: "<<zmAB[0]<<endl;
//                    cout << "IN: "<<t_track.al[4]<<endl;
                      EvaluateClusterPosition_Calorimeter( (int)fmax((float)iclx,(float)icly), sensor, xmABar, ymABar, zmAB, t_track.al[4]); 
//                    cout << ip << "  SENSOR "<<sensor<<endl;
                  }
                  if(isXY){
                      t_track.SetXY(ip,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
                      t_track.SetZ(ip,zmAB[0]);
                  }
                  if(isX ){
                      t_track.SetX( ip,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
                      t_track.SetZ(ip,0.5*(zmAB[1]+zmAB[2]));
                  }
                  if(isY ){
                      t_track.SetY( ip,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
                      t_track.SetZ(ip,0.5*(zmAB[1]+zmAB[2]));
                  }

                  t_track.SetXGood(ip,(sensor == -1 ? 0 : iclx+1),ladder,sensor);
                  t_track.SetYGood(ip,(sensor == -1 ? 0 : icly+1),ladder,sensor);
//                t_track.SetXGood(ip,iclx+1,ladder,sensor);
//                t_track.SetYGood(ip,icly+1,ladder,sensor);
        
                  //=========================================
                  // evaluate other quantities
                  //=========================================
                  float smip,multmax;
                  if(t_track.XGood(ip)){
                      if(isTrk){
                          TrkCluster *cluster = _trk_l1->GetCluster(iclx);
                          float mip = TrkParams::GetMIP(cluster->GetLadder()-1,cluster->view-1);
                          smip    = cluster->GetSignal()/(mip>0.?mip:1.); 
                          multmax = cluster->GetMultiplicity()*10000+cluster->maxs;
                      }else{
                          ExtHit cluster = _cal_cl[iclx];
                          smip    = cluster.signal;//mip
                          multmax = cluster.mult*10000+lrint(cluster.coordPU);
                      }
                      t_track.dedx_x[ip]   = smip/(factor>0.?factor:1.); 
                      t_track.multmaxx[ip] = multmax;
                  }
                  if(t_track.YGood(ip)){
                      if(isTrk){
                          TrkCluster *cluster = _trk_l1->GetCluster(icly);
                          float mip = TrkParams::GetMIP(cluster->GetLadder()-1,cluster->view-1);
                          smip    = cluster->GetSignal()/(mip>0.?mip:1.); 
                          multmax = cluster->GetMultiplicity()*10000+cluster->maxs;
                      }else{
                          ExtHit cluster = _cal_cl[icly];
                          smip    = cluster.signal;//mip
                          multmax = cluster.mult*10000+lrint(cluster.coordPU);
                      }
                      t_track.dedx_y[ip]   = smip/(factor>0.?factor:1.); 
                      t_track.multmaxy[ip] = multmax;
                  }
              }
              ///////////////////////////////////////////////////////// end additional info
 
              // ///////////////////////////////////////////////////////// RE-FIT START
              t_track.Fit(0.,ifail,0);//evaluate again taking into account the track angle
              if( TMath::IsNaN(t_track.chi2) )continue;
              if( ifail != 0 )continue;
              // ///////////////////////////////////////////////////////// RE-FIT FINISH
              if(_debug){
                  cout << endl;
//                cout << " >> al: ";for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
//                cout << " chi2: "<<setw(10)<< t_track.chi2;
//                //            cout << endl;
                  cout<< " TRK chi2 "<<setw(13)<<t_track.chi2;
                  cout <<" | al: ";
                  for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
                  cout <<" | ";
                  cout << " X: ";
                  //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.XGood(ip) ? "o" : "-");
                  for(int ip=0; ip<t_track.nplanes; ip++){
                      if(t_track.GetClusterX_ID(ip)>=0)
                          cout << setw(2)<<t_track.GetClusterX_ID(ip);
                      else
                          cout << "__";
                  }
                  cout << " | ";
                  cout << " Y: ";
                  //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.YGood(ip) ? "o" : "-");
                  for(int ip=0; ip<t_track.nplanes; ip++){
                      if(t_track.GetClusterY_ID(ip)>=0)
                          cout << setw(2)<<t_track.GetClusterY_ID(ip);
                      else
                          cout << "__";
                  }
                  cout << " S:";for(int ip=0; ip<t_track.nplanes; ip++)cout<< t_track.GetSensor(ip)+1;
                  cout << " >> ADD ";
              }

              _caloTj->DoTrack(t_track.al,t_track.zini);
              for(int ip=0; ip<_caloTj->npoint; ip++){
                  t_track.xGF[ip] = _caloTj->x[ip];
                  t_track.yGF[ip] = _caloTj->y[ip];
              }

              //              cout << "ADD "<<endl;
              trackCandidates[mapIndex].insert(pair<int,ExtTrack>(t_track.GetNX()+t_track.GetNY(), ExtTrack(t_track)));
              if(t_track.GetNX()+t_track.GetNY()>5 && t_track.chi2<chi2Min) chi2Min=t_track.chi2;
              iCand++;
          }//end add-candidate condition
            
      }//end first loop on candidates
      if(_debug)cout <<endl<< "n.candidates "<<trackCandidates[mapIndex].size()<<endl;


      //=================================================================
      // ----------------------------------------------------------------
      // chi2 selection
      // ----------------------------------------------------------------
      //=================================================================
      pair <multimap<int,ExtTrack>::iterator, multimap<int,ExtTrack>::iterator> nxny[NEXTVIEWS];
//  for(int inxny=0; inxny<NEXTVIEWS; inxny++){
      bool NMTI =false;//not-minimal track inserted
      for(int inxny=NEXTVIEWS-1; inxny>=nClFixX+nClFixY; inxny--){
          nxny[inxny] = trackCandidates[mapIndex].equal_range(inxny);
//    cout << " --- n "<<inxny<<endl;
          for ( multimap<int,ExtTrack>::iterator it = nxny[inxny].first; it!=nxny[inxny].second; ++it){
              ExtTrack &trackC = (*it).second;
//      cout << " NX+NY "<< trackC.GetNX()+trackC.GetNY() <<" chi2 "<<trackC.chi2<< " R(GV) "<<trackC.GetRigidity()<<endl;
              //            if( trackC.chi2==0 )continue;
              if( (trackC.GetNX()+trackC.GetNY() > 5 && trackC.chi2>chi2Min  ) ||
                  (trackC.GetNX()+trackC.GetNY() == 5 && NMTI ) ||
                  (trackC.GetNX()+trackC.GetNY() == 5 && trackC.chi2>1. ) ||
                  false)continue;      
              
              trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(inxny, ExtTrack(trackC)));
//      cout << " insert "<<endl;
              
              if( trackC.GetNX()+trackC.GetNY() > 5 )NMTI=true;
          }
      }
      mapIndex++;
      
      //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
      if(
          trackCandidates[mapIndex].size()==0 ||
          (trackCandidates[mapIndex].size()>1 && nViewCal < _alg_nViewCal ) ||
          false){
          mapIndex++;
          continue;//increment one calo plane
      }
      //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

      if(_debug){     
    
          cout<< " Selected-track SUMMARY:"<<endl;

//        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
          trkBegin = trackCandidates[mapIndex].begin();
          trkEnd   = trackCandidates[mapIndex].end();
          for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
              ExtTrack &trackCand = (*trkIter).second;
              cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
              cout <<" | al: ";
              for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
              cout <<" | ";
              cout << " X: ";
              //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
              for(int ip=0; ip<trackCand.nplanes; ip++){
                  if(trackCand.GetClusterX_ID(ip)>=0)
                      cout << setw(2)<<trackCand.GetClusterX_ID(ip);
                  else
                      cout << "__";
              }
              cout << " | ";
              cout << " Y: ";
              //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
              for(int ip=0; ip<trackCand.nplanes; ip++){
                  if(trackCand.GetClusterY_ID(ip)>=0)
                      cout << setw(2)<<trackCand.GetClusterY_ID(ip);
                  else
                      cout << "__";
              }
              cout << " S:";for(int ip=0; ip<trackCand.nplanes; ip++)cout<< trackCand.GetSensor(ip)+1;
              cout <<"  MDR="<<1./sqrt(trackCand.coval[4][4])<<endl;
              cout << endl;
              // for(int ip=0; ip<trackCand.nplanes; ip++)
              //        if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";       cout << endl;

          }
      }


      //=================================================================
      // ----------------------------------------------------------------
      // save tracks
      // ----------------------------------------------------------------
      //=================================================================
//      cout << mapIndex << endl;
      if(mapIndex>=maxSize)cout <<" ORRORE!!!! mapIndex"<<mapIndex<<" >= "<<maxSize<<endl;
//      if(mapIndex>0)trackCandidates[mapIndex-1].clear();
      trkBegin = trackCandidates[mapIndex].begin();
      trkEnd   = trackCandidates[mapIndex].end();
      iCand = 0;
      // TClonesArray &tracks = *_trkArray;
      // tracks.Clear();
      for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){        
          ExtTrack trackCand = ((*trkIter).second);
//        cout << iCand << " "<< trackCand.al[4] <<endl;
          //    trackCand.seqno = iCand;
//        trackCand.Dump();
          new(tracks[iCand]) ExtTrack(trackCand);
          iCand++;
      };
      break;


  }// end loop over n.calo planes to be included


}


/**
 * 
 */
void ExtTrkingAlg::ProcessEvent(Bool_t force){
    
    if(_whichAlg == 0)                           ProcessEvent0(force);
    else if(_whichAlg >= 100 && _whichAlg < 144) ProcessEvent1(force);
    else if(_whichAlg >= 200 && _whichAlg < 244) ProcessEvent2(force);
    else{
        cout << " Algorythm id "<< _whichAlg <<" not implemented "<<endl;
    }

};
/**
 * 
 */
void ExtTrkingAlg::Dump(){

  cout << "|--|--|--|--|--|--|--|--|--|--|--|"<<endl;
  cout << "ExtTrkingAlg::Dump() "<<endl<<endl;
  cout << "  algorythm "<<setw(10)<<_whichAlg <<endl;
  cout << "  debug     "<<setw(10)<< _debug <<endl;
  cout << "  TrkLevel1 "<<setw(10)<<_trk_l1 <<endl;
  cout << "  TrkLevel2 "<<setw(10)<<_trk_l2 <<endl;
  cout << "  CaloLevel1"<<setw(10)<<_cal_l1 <<endl;
  cout << "  CaloLevel2"<<setw(10)<<_cal_l2 <<endl;
  cout << "  ToFLevel2 "<<setw(10)<<_tof_l2 <<endl;
  cout << "Pre-selection parameters "<<endl;
  cout << "  nClstrMAX  "<<setw(10)<<_sel_nClstrMAX <<endl;
  cout << "  nPlaneXMIN "<<setw(10)<<_sel_nPlaneXMIN <<endl;
  cout << "  nPlaneYMIN "<<setw(10)<<_sel_nPlaneYMIN <<endl;
  cout << "Algorythm parameters "<<endl;
  cout << "  nClFixX    "<<setw(10)<<_alg_nClFixX <<endl;
  cout << "  nClFixY    "<<setw(10)<<_alg_nClFixY <<endl;
  cout << "  nTrackMAX  "<<setw(10)<<_alg_nTrackMAX <<endl;
  cout << "  nViewCal   "<<setw(10)<<_alg_nViewCal <<endl;
  cout << "|--|--|--|--|--|--|--|--|--|--|--|"<<endl;
  
}


////////////////////////////////////////////////////////////////////////
//
//
//
//
////////////////////////////////////////////////////////////////////////
//  int GetSensor(int view, int sis){
//    if(sis<0||sis>8)return -1;
//    if(view<0||view>1)return -1;
//    if(!view)return (int)(sis/3); 
//    else return sis%3; 
//  };
//  int GetLadder(int view, int sis){ 
//    if(sis<0||sis>8)return -1;
//    if(view<0||view>1)return -1;
//    if(!view)return sis%3;  
//    else return (int)(sis/3); 
//  };
//  int GetSiSensor(int view, int l, int s){
//    if(view<0||view>1)return -1;
//    if(s<0||s>2)return -1;
//    if(!view) return 3*s+l;
//    else      return 3*l+s;
//   };
 
  int CaloStripRoto::GetSensor(){ 
      return GetCaloSensor(GetView(),GetSiSensor()); 
  };
  int CaloStripRoto::GetLadder(){ 
      return GetCaloLadder(GetView(),GetSiSensor()); 
  };

 void CaloStripRoto::ResetAligParams(){
   fPitch = 1.;
   shift[0]=0.; //L 
   shift[1]=0.; //S  
   shift[2]=0.; //Z
   alpha = 0.;
   
//   siSensor = -1;;

 };
/*1
 *
 
 *
 */
bool CaloStripRoto::SensorContains(float x, float y){

    float toll = 0.2;//cm

//    cout << "S -- SetStrip("<<GetLadder()*32+0<<")"<<endl;
    
//    cout << " extst.SetStrip("<<GetLadder()*32+0<<") "<<endl;
    SetStrip(GetLadder()*32+0);
    float xM0,yM0;
    float d0 = GetDistanceTo(x,y,xM0,yM0);

    if( yM0-toll > fYB && yM0 > fYA )return false;
    if( yM0+toll < fYB && yM0 < fYA )return false;

    if(  xM0-toll > fXB && xM0 > fXA )return false;
    if(  xM0+toll < fXB && xM0 < fXA )return false;

//    cout << "D -- SetStrip("<<GetLadder()*32+31<<")"<<endl;
//    cout << " extst.SetStrip("<<GetLadder()*32+31<<") "<<endl;
    SetStrip(GetLadder()*32+31);
    float xM31,yM31;
    float d31 = GetDistanceTo(x,y,xM31,yM31);


    if( d0+d31 > 31*fPitch*0.244+2.*toll)return false;

    return true;

};

 CaloStripRoto::CaloStripRoto(int view, int plane, int sisensor, Bool_t usemechanicalalignement){
     
     SetView(view);
     SetPlane(plane);
     SetSiSensor(sisensor);

     st = CaloStrip(usemechanicalalignement);
     ResetAligParams();
     if(!usemechanicalalignement){
//       cout <<" SETTA I PARAMETRI DI ALLINEAMNETO  -- da implementare"<<endl;
     }
     
//     Set(view,plane,0,GetSiSensor());
 } ///< Default Constructor.


 //         relative to pamela reference system (convenzione emiliano):
 //
 //        y ^
 //          |
 //          | 6 7 8 sisensor
 //          | 3 4 5   
 //          | 0 1 2 
 //          ------------> x 
 //
 //        relative to silicon sensor           
 //
 //        S ^
 //          |   |
 //        2 |   |  <- strip direction
 // sensor 1 |   | 
 //        0 |...|...
 //          ------------> L
 //            0 1 2 
 //            ladder
 //
 //
void CaloStripRoto::Set(Int_t view, Int_t plane, float strip, int sisensor){


    SetView(view);
    SetPlane(plane);


    float size = 8.05; //(cm) larghezza del sensore + zone morte
   
   
   //----------------------------
   // Z-coordinate
   //----------------------------

    float coordCM_strip;
    if( strip >0. && strip < 95.){
        st.Set( view,plane,(int)strip);// istrip
        coordCM_strip = ( view ? st.GetY() : st.GetX() );
        coordCM_strip *= (1-strip+(int)strip);
        st.Set( view,plane,(int)strip+1);// istrip+1
        coordCM_strip += (strip-(int)strip)*( view ? st.GetY() : st.GetX() );
    }else{
//      cout << "CaloStripRoto::Set( "<<view<<","<<plane<<","<<strip<<","<<sisensor<<")"<<endl;
//      cout << " ---> ERROR: 0 =< strip =< 95 " << endl;
        if(strip>=95.)strip=95.;
        if(strip<=0.)strip=0.;
        st.Set( view,plane,(int)strip);// istrip
        coordCM_strip = ( view ? st.GetY() : st.GetX() );       
    }
    float coordCm_Z0 = st.GetZ();
    //----------------------------
    // strip-segment coordinates
    //----------------------------
    float coordCm_LA = coordCM_strip;
    float coordCm_LB = coordCM_strip;
    float coordCm_SA = (-1)*size*1.5; //full calorimeter size
    float coordCm_SB = (+1)*size*1.5; //full calorimeter size
    float coordCm_ZA = coordCm_Z0; 
    float coordCm_ZB = coordCm_Z0; 

  //----------------------------
  // rototraslation
  //----------------------------
   //if sisensor is assigned (0,...8) the strip can be rototraslated
   bool ROTO = (sisensor>=0&&sisensor<9);

   if(ROTO){


       if(         
           GetCaloSiSensor(view,(int)(strip/32),0) != sisensor &&
           GetCaloSiSensor(view,(int)(strip/32),1) != sisensor &&
           GetCaloSiSensor(view,(int)(strip/32),2) != sisensor &&
           true){
           cout << " Strip "<<strip<<" does not belong to sensor "<<sisensor<<endl;
           cout << " input data are inconsistent --- output might be wrong "<<endl;     
       } ;


     int sensor = GetCaloSensor(view, sisensor);//0-1-2
     int ladder = GetCaloLadder(view, sisensor);//0-1-2
     if( ladder != (int)(strip / 32) ){
       cout << " void CaloStripRoto::Set(..) --- ERROR --- strip "<<strip<<" does not match Si-sensor "<<sisensor<<endl;
       return;
     }

     //-----------------------------
     // sensor origin coordinates
     //-----------------------------
     
     st.Set(         view, plane, ladder*32); //strip = 0-32-64
     float coordCm_L0 = ( view ? st.GetY() : st.GetX());  
     st.Set( (int)(!view), plane, sensor*32 );//strip = 0-32-64
     float coordCm_S0 = ( view ? st.GetX() : st.GetY());

     //set sensor size
     coordCm_SA = (-1)*size*0.5 + (sensor-1)*size; 
     coordCm_SB = (+1)*size*0.5 + (sensor-1)*size; 

     //coordinate relative al sensore
     coordCm_LA = coordCm_LA - coordCm_L0; 
     coordCm_SA = coordCm_SA - coordCm_S0;
     coordCm_ZA = coordCm_ZA - coordCm_Z0;
     coordCm_LB = coordCm_LB - coordCm_L0; 
     coordCm_SB = coordCm_SB - coordCm_S0;
     coordCm_ZB = coordCm_ZB - coordCm_Z0;

     //variazione di pitch
     coordCm_LA = coordCm_LA * fPitch;
     coordCm_LB = coordCm_LB * fPitch;

     //rotazione
     float LA = coordCm_LA - alpha * coordCm_SA;
     float SA = coordCm_SA + alpha * coordCm_LA;
     float LB = coordCm_LB - alpha * coordCm_SB;
     float SB = coordCm_SB + alpha * coordCm_LB;

     //traslazione
     coordCm_LA = LA         + (!view ? shift[0] : shift[1]); 
     coordCm_SA = SA         + ( view ? shift[1] : shift[0]); 
     coordCm_ZA = coordCm_ZA + shift[2];
     coordCm_LB = LB         + shift[0];
     coordCm_SB = SB         + shift[1];
     coordCm_ZB = coordCm_ZB + shift[2];

     //back to sistema di riferimento generale
    
     coordCm_LA = coordCm_LA + coordCm_L0; 
     coordCm_SA = coordCm_SA + coordCm_S0;
     coordCm_ZA = coordCm_ZA + coordCm_Z0;
     coordCm_LB = coordCm_LB + coordCm_L0; 
     coordCm_SB = coordCm_SB + coordCm_S0;
     coordCm_ZB = coordCm_ZB + coordCm_Z0;

   }
   
   fXA = (!view ? coordCm_LA : coordCm_SA);
   fYA = ( view ? coordCm_LA : coordCm_SA);
   fZA = coordCm_ZA;

   fXB = (!view ? coordCm_LB : coordCm_SB);
   fYB = ( view ? coordCm_LB : coordCm_SB);
   fZB = coordCm_ZB;


   st.Set(view,plane,(int)strip);// istrip

 };
 //
 //
 //
 //
 //
 //
 float CaloStripRoto::GetSpatialResolution(float def, float degx, float degy,  float beta){
   
   float res_dig = 0.24/sqrt(12.);

   if(def==0.)return res_dig;
   
   float res_ms = 0.;
   
   // Int_t GetView(){return (fView-1);} ///< Get strip view [0-1]
   // Int_t GetPlane(){return (fPlane-1);} ///< Get strip plane [0-21]
   // Int_t GetStrip(){return (fStrip-1);} ///< Get strip number [0-95]
   
   
   float rig = 1/def; if(rig<0)rig = -rig;
   float factor = 0.;
   factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*degx/180.),2);
   factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*degy/180.),2);
   factor += 1.;
   factor = TMath::Sqrt(factor);
   
   int view = st.GetPlane()*2 + 1 - st.GetView()  ;//cluster.view; //0... 43
   
   //  int nW = (int)((view + 1)/2);
   float dW = 0.74*factor;//X0
   float dW_cm = 0.26*factor;//cm
   float dSi_cm = 0.; //cm
   if(view && st.GetPlane()>0) dSi_cm = ( st.GetPlane()%2 ? 0.228 : 0.428);
   dSi_cm *= factor;
   // multiple scattering angle
   float theta = 0.;//ms
   theta = 0.0136/rig;
   theta = theta * sqrt(dW) * (1.+0.038*log(dW));
   float Vt = theta*theta;            //cov(theta,theta)
   float Vy = Vt * dW_cm * dW_cm /3.; //cov(y,y)
   float Vty = 0.87 * Vt * Vy;        //cov(theta,y)
   float Vc = Vy + dSi_cm * dSi_cm * Vt + 2. * dSi_cm * Vty;
   float rms = sqrt( Vc );        
   // if(isY)t_track.resy[ip] = sqrt(rms*rms + t_track.resy[ip]*t_track.resy[ip]);
   // else   t_track.resx[ip] = sqrt(rms*rms + t_track.resx[ip]*t_track.resx[ip]);
   // if(isY)          VMS[nW] = rms*rms;
   // if(_debug)cout <<endl<<view<<" nW "<<nW<<" res(MS) "<<rms<<" factor "<<factor;
   
   if(view>=3){
     cout << " --> NB resolution not accurate for layer "<<view<<endl;;
     cout <<"st.GetPlane() "<<st.GetPlane()<<endl;
     cout <<"st.GetView() "<<st.GetView()<<endl;
   }
   res_ms = rms;
   
   return sqrt(res_ms*res_ms + res_dig*res_dig);
     
   };
/**
 *
 *
 *
 *
 */
float CaloStripRoto::GetDistanceTo(float xP, float yP, float& xM, float &yM ){


    bool isX = fabs(fXA-fXB) < fabs(fYA-fYB);


    if(isX){
        float beta = (fXB-fXA)/(fYB-fYA);
        float alpha = fXA - beta * fYA;
        yM = ( yP + beta*xP - beta*alpha )/(1+beta*beta);
        xM = alpha + beta * yM;
    }else{
        float beta = (fYB-fYA)/(fXB-fXA);
        float alpha = fYA - beta * fXA;
        xM = ( xP + beta*yP - beta*alpha )/(1+beta*beta);
        yM = alpha + beta * xM;
    }
    
    return sqrt((xM-xP)*(xM-xP) + (yM-yP)*(yM-yP));

};