/**
 * \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(true);
   
  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);

    if(_debug)cout << " --->  "<<filep.Data()<<endl;
    std::fstream fs;
    fs.open (filep.Data(), std::fstream::in );

    if(!fs.is_open()){
	
	cout <<" WARNING!!!!!!!!!!!! error opening file "<<filep.Data()<<endl;
	cout <<" No alignment parameters loaded "<<endl;
    }
//    do{
    while(fs.is_open() && !fs.eof()){

	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_l1)return false;
  TrkLevel1 *trkl1 = _trk_l1;//event->GetTrkLevel1();
  int nClstrMAX  = _sel_nClstrMAX; //maximum number of cluster
  /////////////////////////////////////////////////////////////////////
  /// dump selection
  /////////////////////////////////////////////////////////////////////
  if(_debug){
    cout << " n.clusters "<<trkl1->nclstr()<<endl;      
  }  
  if(trkl1->nclstr() > nClstrMAX) return false; 
  /////////////////////////////////////////////////////////////////////
  /// count number of hit planes
  /////////////////////////////////////////////////////////////////////
  int nPlaneXMIN = _sel_nPlaneXMIN; 
  int nPlaneYMIN = _sel_nPlaneYMIN;
  int nHitX[] = {0,0,0,0,0,0};
  int nHitY[] = {0,0,0,0,0,0};
  for(int ic=0; ic<trkl1->nclstr(); ic++){      
      if      (trkl1->GetCluster(ic)->IsX() ) nHitX[trkl1->GetCluster(ic)->GetPlane()-1]++;    
      else if (trkl1->GetCluster(ic)->IsY() ) nHitY[trkl1->GetCluster(ic)->GetPlane()-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.x-hit planes  "<<nPlaneX<<endl;
      cout << " n.y-hit planes  "<<nPlaneY<<endl;
  }  
  if( nPlaneX<nPlaneXMIN ) return false; 
  if( nPlaneY<nPlaneYMIN ) return false; 


//   if(!_trk_l2)return true;
//   TrkLevel2 *trkl2 = _trk_l2;//event->GetTrkLevel2();
//   /////////////////////////////////////////////////////////////////////
//   /// 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.y-singles "<<trkl2->nclsy()<<endl;
//   }
//   if(trkl2->nclsy()+trkl2->nclsx() > nClstrMAX) return false; 
//   /////////////////////////////////////////////////////////////////////
//   /// count number of hit planes
//   /////////////////////////////////////////////////////////////////////
//   int nPlaneXMIN = _sel_nPlaneXMIN; 
//   int nPlaneYMIN = _sel_nPlaneYMIN;
//   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.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; 

  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));

};