/[PAMELA software]/DarthVader/TrackerLevel2/src/ExtTrkingAlg.cpp

Diff of /DarthVader/TrackerLevel2/src/ExtTrkingAlg.cpp

Parent Directory | Revision Log | View Patch Patch

-revision 1.1 by pam-fi,
Thu Feb 27 11:24:43 2014 UTC
+revision 1.5 by pam-ts,
Thu Jul 24 12:39:50 2014 UTC
 Line 9
  using namespace std;
  /**
   * Constructor
+  *
+  * @param id Algorythm identifier
+  *
+  * id = 0 The algorythm looks for tracks with >= 3x2y hits on the tracker plane.
+  *        It fills a TClonesArray of TrkTrack objects.
+  * id = 1 The first plane (y) of the calorimeter is also used.
+  *        It fills a TClonesArray of ExtTrack objects.
+  * id = 2 The first two planes (y+x) of the calorimeter are also used.
+  *        It fills a TClonesArray of ExtTrack objects.
+  *
   */
  ExtTrkingAlg::ExtTrkingAlg(Int_t id){
    _whichAlg = id;
-   SetDebug(false);
-   double pars[] = {20.,3.,2.};
+   SetDebug(true);
-   SetSelectionParams(pars);
-   double para[] = {3.,2.,400.};
-   SetAlgorythmParams(para);
    if(id == 0){
+ //    if(_debug){
-     cout << "ExtTrkingAlg::ExtTrkingAlg("<<id<<")"<<endl;
+       cout << "ExtTrkingAlg::ExtTrkingAlg("<<id<<")"<<endl;
-     cout << "Creating array of TrkTrack objects "<<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{
+   // }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
-     cout << "ExtTrkingAlg(Int_t id) -- algorythm id="<<id<<" not valid "<<endl;
+     if((int)(id/100)==1)_alg_nViewCal = id - 100;  // algorythm parameter: n.calorimeter planes included
-     cout << "--> Track array not created "<<endl;
+     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();
  };
-Line 43 
 ExtTrkingAlg::ExtTrkingAlg(Int_t id){
+Line 222 
 ExtTrkingAlg::ExtTrkingAlg(Int_t id){
  /**
   * Reset the algorythm flags and output
   */
- void ExtTrkingAlg::Reset(){
+ // void ExtTrkingAlg::Reset(){
-   if(_trkArray)_trkArray->Clear("C");
+ //   if(_trkArray)_trkArray->Clear("C");
- };
+ // };
  /**
   * Clear the event and reset the algorythm
-Line 64 
 void ExtTrkingAlg::Clear(Option_t* optio
+Line 241 
 void ExtTrkingAlg::Clear(Option_t* optio
    SetCaloLevel1();
    SetCaloLevel2();
-   Reset();
+   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
   */
-Line 85 
 void ExtTrkingAlg::SetAlgorythmParams(do
+Line 274 
 void ExtTrkingAlg::SetAlgorythmParams(do
    _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];     //
  };
  /**
-Line 93 
 void ExtTrkingAlg::SetAlgorythmParams(do
+Line 284 
 void ExtTrkingAlg::SetAlgorythmParams(do
   */
  TClonesArray* ExtTrkingAlg::GetTrackArray(Bool_t reset){
-   if( reset ) Reset();
+   if( reset && _trkArray) _trkArray->Clear("C");
    //  if( !_procDone )ProcessEvent();
    return _trkArray;
-Line 104 
 TClonesArray* ExtTrkingAlg::GetTrackArra
+Line 295 
 TClonesArray* ExtTrkingAlg::GetTrackArra
   */
  Bool_t ExtTrkingAlg::CheckEvent(){
-   if(!_trk_l2)return true;
+   if(!_trk_l1)return false;
+   TrkLevel1 *trkl1 = _trk_l1;//event->GetTrkLevel1();
-   TrkLevel2 *trkl2 = _trk_l2;//event->GetTrkLevel2();
    int nClstrMAX  = _sel_nClstrMAX; //maximum number of cluster
-   int nPlaneXMIN = _sel_nPlaneXMIN;
+   /////////////////////////////////////////////////////////////////////
-   int nPlaneYMIN = _sel_nPlaneYMIN;
+   /// 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 ix=0; ix<trkl2->nclsx(); ix++)nHitX[trkl2->GetSingletX(ix)->plane-1]++;
+   for(int ic=0; ic<trkl1->nclstr(); ic++){
-   for(int iy=0; iy<trkl2->nclsy(); iy++)nHitY[trkl2->GetSingletY(iy)->plane-1]++;
+       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
    /////////////////////////////////////////////////////////////////////
-Line 124 
 Bool_t ExtTrkingAlg::CheckEvent(){
+Line 323 
 Bool_t ExtTrkingAlg::CheckEvent(){
    for(int ip=0; ip<6; ip++)if(nHitX[ip])nPlaneX++;
    int nPlaneY=0;
    for(int ip=0; ip<6; ip++)if(nHitY[ip])nPlaneY++;
    /////////////////////////////////////////////////////////////////////
    /// dump selection
    /////////////////////////////////////////////////////////////////////
    if(_debug){
-     cout << " n.tracks "<<trkl2->GetNTracks()<<endl;
+       cout << " n.x-hit planes  "<<nPlaneX<<endl;
-     cout << " n.x-singles "<<trkl2->nclsx()<<endl;
+       cout << " n.y-hit planes  "<<nPlaneY<<endl;
-     cout << " n.y-singles "<<trkl2->nclsy()<<endl;
+   }
-     cout << " n.x-planes with singles "<<nPlaneX<<endl;
+   if( nPlaneX<nPlaneXMIN ) return false;
-     cout << " n.y-planes with singles "<<nPlaneY<<endl;
+   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 && nPlaneX<nPlaneXMIN) return false;
-   if(trkl2->GetNTracks()==0 && nPlaneY<nPlaneYMIN) return false;
+ //   if(trkl2->GetNTracks()==0 && nPlaneY<nPlaneYMIN) return false;
-   if(trkl2->nclsy()+trkl2->nclsx() > nClstrMAX) return false;
    return true;
-Line 151 
 Bool_t ExtTrkingAlg::CheckEvent(){
+Line 386 
 Bool_t ExtTrkingAlg::CheckEvent(){
  /**
   * Process the event
   */
- void ExtTrkingAlg::ProcessEvent(Bool_t force){
+ void ExtTrkingAlg::ProcessEvent0(Bool_t force){
    if(_debug)cout << " |---------------------------------------------------| "<<endl;
-   if(_debug)cout << " void ExtTrkingAlg::ProcessEvent("<<force<<") "<<endl;
+   if(_debug)cout << " void ExtTrkingAlg::ProcessEvent0("<<force<<") "<<endl;
    if(_debug && !_trk_l1)cout << " Missing TrkLevel1 object --- ciao! "<<endl;
-Line 175 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
+Line 410 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
    // tracking
    //===========================
-   TClonesArray &tracks = *_trkArray;
+   // TClonesArray &tracks = *_trkArray;
-   tracks.Clear();
+   // tracks.Clear();
    multimap<int,TrkTrack> trackCandidates[12]; // map: last filled view vs TrkTrack
    int mapIndex = 0;
    multimap<int,TrkTrack>::iterator trkBegin;
-Line 304 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
+Line 541 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
    }//end iterations
-   if( trackCandidates[mapIndex].size() > nTrackMAX ){
+   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
    }
-Line 316 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
+Line 553 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
        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
    // -------------------------------------------------------------
-Line 391 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
+Line 628 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
      // -----------------------------------
      // evaluate calorimeter variables
      // -----------------------------------
-     float caloCoord[2][22][96];
+     // float caloCoord[2][22][96];
-     float caloZ[44];
+     // float caloZ[44];
-     for(int view=0; view<2; view++){
+     // for(int view=0; view<2; view++){
-       for(int plane=0; plane<22; plane++){
+     //   for(int plane=0; plane<22; plane++){
-         for(int strip=0; strip<96; strip++){
+     //  for(int strip=0; strip<96; strip++){
-           CaloStrip st = CaloStrip(false);
+     //    CaloStrip st = CaloStrip(false);
-           st.Set(view,plane,strip);
+     //    st.Set(view,plane,strip);
-           caloCoord[view][plane][strip]=(view ? st.GetY() : st.GetX());
+     //    caloCoord[view][plane][strip]=(view ? st.GetY() : st.GetX());
-           caloZ[plane*2+(view?0:1)]=st.GetZ();
+     //    caloZ[plane*2+(view?0:1)]=st.GetZ();
-         }
+     //  }
-       }
+     //   }
-     }
+     // }
-     Trajectory caloTj = Trajectory(44,caloZ);
+     // Trajectory caloTj = Trajectory(44,caloZ);
      //  for(int iz=0; iz<44;iz++)cout<<caloZ[iz]<<endl;
      //test calo
      //    float caloEvent[2][22][96];
-Line 417 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
+Line 654 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
      for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
        TrkTrack trackCand = ((*trkIter).second);
-       trackCand.DoTrack(&caloTj);//integrate the trajectory through the calorimeter
+       trackCand.DoTrack(_caloTj);//integrate the trajectory through the calorimeter
        float chi2 = 0.;
        int  nhit = 0;
-Line 431 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
+Line 668 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
          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)]);
+         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++;
-Line 748 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
+Line 985 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
    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;
-Line 759 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
+Line 998 
 void ExtTrkingAlg::ProcessEvent(Bool_t f
  };
+ /**
+  * 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::Dump(){
+ void ExtTrkingAlg::FillClusterMap(multimap<int,int> &map,CaloLevel1* l1,Int_t vOffset){
-   cout << "|--|--|--|--|--|--|--|--|--|--|--|"<<endl;
+   if(_debug)cout<<"-------- CALO CLUSTERs -------- on VIEW < "<<_alg_nViewCal<<endl;
-   cout << "ExtTrkingAlg::Dump() "<<endl<<endl;
-   cout << "  algorythm "<<setw(10)<<_whichAlg <<endl;
+   if(_alg_nViewCal<=0)return;
-   cout << "  debug     "<<setw(10)<< _debug <<endl;
-   cout << "  TrkLevel1 "<<setw(10)<<_trk_l1 <<endl;
+   _cal_cl.clear();
-   cout << "  TrkLevel2 "<<setw(10)<<_trk_l2 <<endl;
-   cout << "  CaloLevel1"<<setw(10)<<_cal_l1 <<endl;
+   if(!l1)return;
-   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;
+   // CALORIMETER CLUSTERING
-   cout << "  nPlaneXMIN "<<setw(10)<<_sel_nPlaneXMIN <<endl;
+   // ==========================
-   cout << "  nPlaneYMIN "<<setw(10)<<_sel_nPlaneYMIN <<endl;
+   CaloStrip st = CaloStrip(false);
-   cout << "Algorythm parameters "<<endl;
+   ExtHit *caloCl=NULL;
-   cout << "  nClFixX    "<<setw(10)<<_alg_nClFixX <<endl;
+   Int_t nCl = 0;
-   cout << "  nClFixY    "<<setw(10)<<_alg_nClFixY <<endl;
+   for(Int_t ih=0; ih<l1->istrip; ih++){//loop over calorimeter hits
-   cout << "  nTrackMAX  "<<setw(10)<<_alg_nTrackMAX <<endl;
+     Int_t view = -1;  // 0-x 1-y
-   cout << "|--|--|--|--|--|--|--|--|--|--|--|"<<endl;
+     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));
+ };

 Legend:



Removed from v.1.1
 


changed lines


 
Added in v.1.5
 Legend:



Removed from v.1.1
 


changed lines


 
Added in v.1.5
-Removed from v.1.1
+Added in v.1.5

	ViewVC Help
Powered by ViewVC 1.1.23