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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by pam-fi, Thu Feb 27 11:24:43 2014 UTC revision 1.4 by pam-ts, Mon Jun 23 08:34:18 2014 UTC
# Line 9  Line 9 
9  using namespace std;  using namespace std;
10  /**  /**
11   * Constructor   * Constructor
12     *
13     * @param id Algorythm identifier
14     *
15     * id = 0 The algorythm looks for tracks with >= 3x2y hits on the tracker plane.
16     *        It fills a TClonesArray of TrkTrack objects.
17     * id = 1 The first plane (y) of the calorimeter is also used.
18     *        It fills a TClonesArray of ExtTrack objects.
19     * id = 2 The first two planes (y+x) of the calorimeter are also used.
20     *        It fills a TClonesArray of ExtTrack objects.
21     *
22   */   */
23  ExtTrkingAlg::ExtTrkingAlg(Int_t id){  ExtTrkingAlg::ExtTrkingAlg(Int_t id){
24    
# Line 16  ExtTrkingAlg::ExtTrkingAlg(Int_t id){ Line 26  ExtTrkingAlg::ExtTrkingAlg(Int_t id){
26    
27    SetDebug(false);    SetDebug(false);
28    
   double pars[] = {20.,3.,2.};  
   SetSelectionParams(pars);  
   
   
   double para[] = {3.,2.,400.};  
   SetAlgorythmParams(para);  
29        
30    if(id == 0){    if(id == 0){
31        if(_debug){
32      cout << "ExtTrkingAlg::ExtTrkingAlg("<<id<<")"<<endl;        cout << "ExtTrkingAlg::ExtTrkingAlg("<<id<<")"<<endl;
33      cout << "Creating array of TrkTrack objects "<<endl;        cout << "Creating array of TrkTrack objects "<<endl;
34          cout << "WARNING!!! tracking not accurate!! bug not fixed yet... "<<endl;
35        }
36      _trkArray = new TClonesArray("TrkTrack");      _trkArray = new TClonesArray("TrkTrack");
37        //
38        _sel_nClstrMAX  = 20. ; // selection parameter: maximum number of cluster
39        _sel_nPlaneXMIN = 3.; // selection parameter: minimum number of hit x-views
40        _sel_nPlaneYMIN = 2.; // selection parameter: minimum number of hit y-views
41        //
42        _alg_nClFixX   = 3.;     // algorythm parameter: n.hits required on X view
43        _alg_nClFixY   = 2.;     // algorythm parameter:n.hits required on Y view
44        _alg_nTrackMAX = 400.;     // algorythm parameter: maximum num. of track candidates
45        _alg_nViewCal= 0; // ignored
46        //    _alg_nViewCalFix= 0; // ignored
47    
48        NEXTVIEWS  = 12 +_alg_nViewCal;
49        NEXTPLANES = 6  +_alg_nViewCal;
50        _extTrack = NULL;//ignored
51        _zMech = NULL;//ignored
52        
53    
54    }else{    // }else if (id == 1 ){
55      }else if ( (id >= 100 && id < 144) ||
56                 (id >= 200 && id < 244) ||
57                 false){ //add n calo planes
58        if(_debug){
59          cout << "ExtTrkingAlg::ExtTrkingAlg("<<id<<")"<<endl;
60          cout << "Creating array of ExtTrack objects "<<endl;
61        }
62        _trkArray = new TClonesArray("ExtTrack");
63        //
64        _alg_nClFixX   = 3.;     // algorythm parameter: n.hits required on X view
65        _alg_nClFixY   = 2.;     // algorythm parameter:n.hits required on Y view
66        _alg_nTrackMAX = 1000.;     // algorythm parameter: maximum num. of track candidates
67        // _alg_nViewCal = 1;        // algorythm parameter: n.calorimeter planes included
68        //_alg_nViewCal = 1;        // algorythm parameter: n.calorimeter planes included
69        
70        if((int)(id/100)==1)_alg_nViewCal = id - 100;  // algorythm parameter: n.calorimeter planes included
71        if((int)(id/200)==1)_alg_nViewCal = id - 200;  // algorythm parameter: n.calorimeter planes included
72    
73        //
74        _sel_nClstrMAX  = 15.;//10. ; // selection parameter: maximum number of cluster
75        _sel_nPlaneXMIN = (_alg_nViewCal ? 2 : 3);   // selection parameter: minimum number of hit x-views
76        _sel_nPlaneYMIN = 2 - (int)(_alg_nViewCal+1)/2 ;   // selection parameter: minimum number of hit y-views
77    
78        NEXTVIEWS  = 12 +_alg_nViewCal;
79        NEXTPLANES = 6  +_alg_nViewCal;
80        _extTrack = new ExtTrack(NEXTPLANES);
81        _zMech = new float[NEXTPLANES];
82        
83        // -----------------------------------
84        // set Z-coord of tracker planes
85        // -----------------------------------  
86        _zMech[0]=ZTRK1;
87        _zMech[1]=ZTRK2;
88        _zMech[2]=ZTRK3;
89        _zMech[3]=ZTRK4;
90        _zMech[4]=ZTRK5;
91        _zMech[5]=ZTRK6;
92        // _extTrack->SetZ(0,ZTRK1);
93        // _extTrack->SetZ(1,ZTRK2);
94        // _extTrack->SetZ(2,ZTRK3);
95        // _extTrack->SetZ(3,ZTRK4);
96        // _extTrack->SetZ(4,ZTRK5);
97        // _extTrack->SetZ(5,ZTRK6);
98            
99      cout << "ExtTrkingAlg(Int_t id) -- algorythm id="<<id<<" not valid "<<endl;      _caloStripRoto.clear();
100      cout << "--> Track array not created "<<endl;      for(int iv=0; iv<_alg_nViewCal; iv++){ //loop over calorimater tracking planes
101          for(int is=0; is<9; is++){ //loop over calorimeter sensors
102              CaloStripRoto s(1-iv%2, (int)(iv/2), is, false);//flight 0-order align        
103              s.ResetAligParams();//reset 1-order align parameters  
104              _caloStripRoto.push_back(s) ;
105          }
106        }
107        //------------------------------------------
108        // read alignment parameters
109        //------------------------------------------
110        if(_debug) cout << "Reading calorimeter alignment parameters"<<endl;
111        const char *pamca = gSystem->Getenv("PAM_CALIB");
112        TString filep = "/trk-param/align_param_calo-0/CaloAlignParams.txt";
113        filep.Prepend(pamca);
114        std::fstream fs;
115        fs.open (filep.Data(), std::fstream::in );
116        do{
117            int view,plane,sensor;
118            float par[5],err[5];
119            fs>>view >> plane>> sensor;
120            for(int i=0; i<5; i++)fs>>par[i]>>err[i];
121            
122            
123            if(plane*2+1-view<_alg_nViewCal){
124    
125    
126                int index = (plane*2+1-view)*9 + sensor;
127                if(_debug){
128                  cout<<index<<" >> ";
129                  cout<<view << " "<<plane<<" "<<sensor<<"          ";
130                  for(int i=0; i<5; i++)cout<<setw(10)<<par[i]<<" "<<setw(10)<<err[i]<<" ";
131                  cout<<endl;
132                }
133    
134                _caloStripRoto[index].SetAligParams__fPitch(par[0]);
135                _caloStripRoto[index].SetAligParams__shift(par+1);
136                _caloStripRoto[index].SetAligParams__alpha(par[4]);
137    //          _caloStripRoto[index].SetSiSensor(sensor);
138    //          _caloStripRoto[index].SetView(view);
139    //          _caloStripRoto[index].SetPlane(plane);
140    
141            }
142        }while(!fs.eof());
143        fs.close();
144    
145        // -----------------------------------
146        // set Z-coord of calorimeter planes
147        // -----------------------------------      
148        for(int plane=0; plane<22; plane++){
149          for(int view=1; view>=0; view--){
150            int strip = 0;
151            CaloStrip st = CaloStrip(false);//flight align
152            st.Set(view,plane,strip);      
153            int index = 2*plane + !view;
154            if(index >= _alg_nViewCal)break;
155            // _extTrack->SetZ(6+index,st.GetZ());
156            _zMech[6+index]=st.GetZ();
157          }
158        }
159    
160        for(int ip =0 ; ip<_extTrack->nplanes; ip++ )_extTrack->SetZ(ip,_zMech[ip]);
161    
162        if ( _debug ){
163          cout <<" Extended-track Z-coordinates: "<<endl;
164          for(int ip =0 ; ip<_extTrack->nplanes; ip++ )cout <<  _extTrack->zm[ip]<<endl;
165        }
166      }else{
167        if ( _debug ){
168          cout << "ExtTrkingAlg(Int_t id) -- algorythm id="<<id<<" not valid "<<endl;
169          cout << "--> Track array not created "<<endl;
170        }
171        _sel_nClstrMAX  = 0.;
172        _sel_nPlaneXMIN = 0.;  
173        _sel_nPlaneYMIN = 0.;  
174        //
175        _alg_nClFixX   = 0.;  
176        _alg_nClFixY   = 0.;  
177        _alg_nTrackMAX = 0.;  
178        _alg_nViewCal = 0;    
179        //    _alg_nViewCalFix= 0;
180    
181        NEXTVIEWS  = 12 +_alg_nViewCal;
182        NEXTPLANES = 6  +_alg_nViewCal;
183        _extTrack = NULL;//ignored
184    
185      }
186    
187      CaloStrip st = CaloStrip(false);
188      for(int view=0; view<2; view++){
189        for(int plane=0; plane<22; plane++){
190          for(int strip=0; strip<96; strip++){
191            st.Set(view,plane,strip);
192            _caloCoord[view][plane][strip]=(view ? st.GetY() : st.GetX());
193            _caloZ[plane*2+(view?0:1)]=st.GetZ();
194          }
195        }
196    }    }
197      _caloTj = new Trajectory(44,_caloZ);
198    //  _caloTj->Dump();
199    
200      int    ngf = TrkParams::nGF;
201      float *zgf = TrkParams::zGF;
202      _tgf = new Trajectory(ngf,zgf);
203    //  _tgf->Dump();
204    
205    Clear();    Clear();
206        
207  };  };
# Line 43  ExtTrkingAlg::ExtTrkingAlg(Int_t id){ Line 209  ExtTrkingAlg::ExtTrkingAlg(Int_t id){
209  /**  /**
210   * Reset the algorythm flags and output   * Reset the algorythm flags and output
211   */   */
212  void ExtTrkingAlg::Reset(){  // void ExtTrkingAlg::Reset(){
   
   if(_trkArray)_trkArray->Clear("C");  
213    
214      //   if(_trkArray)_trkArray->Clear("C");
215    
216  };  // };
217    
218  /**  /**
219   * Clear the event and reset the algorythm   * Clear the event and reset the algorythm
# Line 64  void ExtTrkingAlg::Clear(Option_t* optio Line 228  void ExtTrkingAlg::Clear(Option_t* optio
228    SetCaloLevel1();    SetCaloLevel1();
229    SetCaloLevel2();    SetCaloLevel2();
230    
231    Reset();    if(_trkArray)_trkArray->Clear("C");
232    
233    };  
234    /**
235     * DElete all
236     */
237    
238    void ExtTrkingAlg::Delete(){
239      if(_trkArray) delete _trkArray;
240      if(_caloTj) delete _caloTj;
241      if(_tgf) delete _tgf;
242      if(_zMech)delete [] _zMech;
243    }
244    
245    
 };  
246  /**  /**
247   * Set selection parameters   * Set selection parameters
248   */   */
# Line 85  void ExtTrkingAlg::SetAlgorythmParams(do Line 261  void ExtTrkingAlg::SetAlgorythmParams(do
261    _alg_nClFixX   = (Int_t)par[0];     // algorythm parameter: n.hits required on X view    _alg_nClFixX   = (Int_t)par[0];     // algorythm parameter: n.hits required on X view
262    _alg_nClFixY   = (Int_t)par[1];     // algorythm parameter:n.hits required on X view    _alg_nClFixY   = (Int_t)par[1];     // algorythm parameter:n.hits required on X view
263    _alg_nTrackMAX = (Int_t)par[2];     // algorythm parameter: maximum num. of track candidates    _alg_nTrackMAX = (Int_t)par[2];     // algorythm parameter: maximum num. of track candidates
264      _alg_nViewCal=  (Int_t)par[3];     // algorythm parameter: n. calorimeter planes included
265      //  _alg_nViewCalFix=  (Int_t)par[4];     //
266    
267  };  };
268    
269  /**  /**
# Line 93  void ExtTrkingAlg::SetAlgorythmParams(do Line 271  void ExtTrkingAlg::SetAlgorythmParams(do
271   */   */
272  TClonesArray* ExtTrkingAlg::GetTrackArray(Bool_t reset){  TClonesArray* ExtTrkingAlg::GetTrackArray(Bool_t reset){
273        
274    if( reset ) Reset();    if( reset && _trkArray) _trkArray->Clear("C");
275    //  if( !_procDone )ProcessEvent();    //  if( !_procDone )ProcessEvent();
276    
277    return _trkArray;    return _trkArray;
# Line 151  Bool_t ExtTrkingAlg::CheckEvent(){ Line 329  Bool_t ExtTrkingAlg::CheckEvent(){
329  /**  /**
330   * Process the event   * Process the event
331   */   */
332  void ExtTrkingAlg::ProcessEvent(Bool_t force){  void ExtTrkingAlg::ProcessEvent0(Bool_t force){
333    
334    if(_debug)cout << " |---------------------------------------------------| "<<endl;    if(_debug)cout << " |---------------------------------------------------| "<<endl;
335    if(_debug)cout << " void ExtTrkingAlg::ProcessEvent("<<force<<") "<<endl;    if(_debug)cout << " void ExtTrkingAlg::ProcessEvent0("<<force<<") "<<endl;
336    
337        
338    if(_debug && !_trk_l1)cout << " Missing TrkLevel1 object --- ciao! "<<endl;    if(_debug && !_trk_l1)cout << " Missing TrkLevel1 object --- ciao! "<<endl;
# Line 175  void ExtTrkingAlg::ProcessEvent(Bool_t f Line 353  void ExtTrkingAlg::ProcessEvent(Bool_t f
353    // tracking    // tracking
354    //===========================    //===========================
355    
356    TClonesArray &tracks = *_trkArray;    // TClonesArray &tracks = *_trkArray;
357    tracks.Clear();    // tracks.Clear();
358    
359    
360    multimap<int,TrkTrack> trackCandidates[12]; // map: last filled view vs TrkTrack    multimap<int,TrkTrack> trackCandidates[12]; // map: last filled view vs TrkTrack
361    int mapIndex = 0;    int mapIndex = 0;
362    multimap<int,TrkTrack>::iterator trkBegin;    multimap<int,TrkTrack>::iterator trkBegin;
# Line 304  void ExtTrkingAlg::ProcessEvent(Bool_t f Line 484  void ExtTrkingAlg::ProcessEvent(Bool_t f
484    }//end iterations    }//end iterations
485        
486    
487    if( trackCandidates[mapIndex].size() > nTrackMAX ){    if( trackCandidates[mapIndex].size() > (uint)nTrackMAX ){ // EM,  warning: comparison between signed and unsigned integer expressions [-Wsign-compare]
488      if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<" > "<<nTrackMAX<<endl;      if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<" > "<<nTrackMAX<<endl;
489      return;//to many candidates      return;//to many candidates
490    }    }
# Line 316  void ExtTrkingAlg::ProcessEvent(Bool_t f Line 496  void ExtTrkingAlg::ProcessEvent(Bool_t f
496        cout << " >> ";        cout << " >> ";
497        cout << " X: ";        cout << " X: ";
498        for(int ip=0; ip<6; ip++)        for(int ip=0; ip<6; ip++)
499          if(trackCand.GetClusterX_ID(ip)>=0)cout << trackCand.GetClusterX_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";          if(trackCand.GetClusterX_ID(ip)>=0)cout << trackCand.GetClusterX_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";
500        cout << " Y: ";        cout << " Y: ";
501        for(int ip=0; ip<6; ip++)        for(int ip=0; ip<6; ip++)
502          if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";              if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";    
503        cout << endl;        cout << endl;
504      }      }
505    }    }
506    if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;    if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
507    
508  //==============================================================    //==============================================================
509    // -------------------------------------------------------------    // -------------------------------------------------------------
510    // first selection of track candidates, based on chi2    // first selection of track candidates, based on chi2
511    // -------------------------------------------------------------    // -------------------------------------------------------------
# Line 391  void ExtTrkingAlg::ProcessEvent(Bool_t f Line 571  void ExtTrkingAlg::ProcessEvent(Bool_t f
571      // -----------------------------------      // -----------------------------------
572      // evaluate calorimeter variables      // evaluate calorimeter variables
573      // -----------------------------------      // -----------------------------------
574      float caloCoord[2][22][96];      // float caloCoord[2][22][96];
575      float caloZ[44];      // float caloZ[44];
576      for(int view=0; view<2; view++){      // for(int view=0; view<2; view++){
577        for(int plane=0; plane<22; plane++){      //   for(int plane=0; plane<22; plane++){
578          for(int strip=0; strip<96; strip++){      //  for(int strip=0; strip<96; strip++){
579            CaloStrip st = CaloStrip(false);      //    CaloStrip st = CaloStrip(false);
580            st.Set(view,plane,strip);      //    st.Set(view,plane,strip);
581            caloCoord[view][plane][strip]=(view ? st.GetY() : st.GetX());      //    caloCoord[view][plane][strip]=(view ? st.GetY() : st.GetX());
582            caloZ[plane*2+(view?0:1)]=st.GetZ();      //    caloZ[plane*2+(view?0:1)]=st.GetZ();
583          }      //  }
584        }      //   }
585      }      // }
586      Trajectory caloTj = Trajectory(44,caloZ);      // Trajectory caloTj = Trajectory(44,caloZ);
587      //  for(int iz=0; iz<44;iz++)cout<<caloZ[iz]<<endl;      //  for(int iz=0; iz<44;iz++)cout<<caloZ[iz]<<endl;
588      //test calo      //test calo
589      //    float caloEvent[2][22][96];      //    float caloEvent[2][22][96];
# Line 417  void ExtTrkingAlg::ProcessEvent(Bool_t f Line 597  void ExtTrkingAlg::ProcessEvent(Bool_t f
597      for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){              for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){        
598        TrkTrack trackCand = ((*trkIter).second);          TrkTrack trackCand = ((*trkIter).second);  
599                            
600        trackCand.DoTrack(&caloTj);//integrate the trajectory through the calorimeter        trackCand.DoTrack(_caloTj);//integrate the trajectory through the calorimeter
601                            
602        float chi2 = 0.;        float chi2 = 0.;
603        int  nhit = 0;        int  nhit = 0;
# Line 431  void ExtTrkingAlg::ProcessEvent(Bool_t f Line 611  void ExtTrkingAlg::ProcessEvent(Bool_t f
611          if(plane<0)continue;          if(plane<0)continue;
612          if(strip<0)continue;          if(strip<0)continue;
613          if(mip<=0)continue;          if(mip<=0)continue;
614          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)]);
615          if(plane<11) chi2+= mip * TMath::Power(dxy,2.);          if(plane<11) chi2+= mip * TMath::Power(dxy,2.);
616          //            caloZ[plane*2+(view?0:1)]=st.GetZ();          //            caloZ[plane*2+(view?0:1)]=st.GetZ();
617          if( fabs(dxy) < 2. && plane<11)nhit++;          if( fabs(dxy) < 2. && plane<11)nhit++;
# Line 748  void ExtTrkingAlg::ProcessEvent(Bool_t f Line 928  void ExtTrkingAlg::ProcessEvent(Bool_t f
928    trkBegin = trackCandidates[mapIndex].begin();    trkBegin = trackCandidates[mapIndex].begin();
929    trkEnd   = trackCandidates[mapIndex].end();    trkEnd   = trackCandidates[mapIndex].end();
930    iCand = 0;    iCand = 0;
931      TClonesArray &tracks = *_trkArray;
932      tracks.Clear();
933    for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){                for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){            
934      TrkTrack trackCand = ((*trkIter).second);      TrkTrack trackCand = ((*trkIter).second);
935      trackCand.seqno = iCand;      trackCand.seqno = iCand;
# Line 759  void ExtTrkingAlg::ProcessEvent(Bool_t f Line 941  void ExtTrkingAlg::ProcessEvent(Bool_t f
941    
942    
943  };  };
944    /**
945     * Fill cluster map with tracker clusters
946     */
947    void ExtTrkingAlg::FillClusterMap(multimap<int,int> &map,TrkLevel1* l1,Int_t vOffset){
948    
949      if(_debug)cout<<"-------- TRK CLUSTERs --------"<<endl;
950    
951      // _trk_cl.clear();
952    
953      if(!l1)return;
954    
955      ExtHit *trkCl = new ExtHit();
956    
957      for(int icl=0; icl<l1->nclstr(); icl++){
958    
959        // trkCl->Reset();
960    
961        // l1->GetCluster(icl)->GetLevel1Struct();// to use F77 routines //PERICOLOSO
962        // trkCl->coordPU = (float)(l1->GetCluster(icl)->maxs) - 1. + l1->GetCluster(icl)->GetCOG();
963        // trkCl->coordCm = 1000.; //TEMPORANEO forse non serve
964        // trkCl->resCm = 1000.;//TEMPORANEO
965        // trkCl->signal = l1->GetCluster(icl)->GetSignal();
966        // trkCl->mult = l1->GetCluster(icl)->GetMultiplicity();
967        // trkCl->view = l1->GetCluster(icl)->view-1;
968    
969        // _trk_cl.push_back(ExtHit(*trkCl));
970    
971        map.insert( pair<int,int>( l1->GetCluster(icl)->view-1 + vOffset, icl ));
972        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;
973      }
974      
975      // l1->GetLevel1Struct(); //necessario per ritornare al contenuto originale del common
976    
977      if(trkCl)delete trkCl;
978    };
979    
980  /**  /**
981   *   * Fill cluster map with calorimeter clusters
982   */   */
983  void ExtTrkingAlg::Dump(){  void ExtTrkingAlg::FillClusterMap(multimap<int,int> &map,CaloLevel1* l1,Int_t vOffset){
984    
985      if(_debug)cout<<"-------- CALO CLUSTERs -------- on VIEW < "<<_alg_nViewCal<<endl;
986    
987      if(_alg_nViewCal<=0)return;
988    
989      _cal_cl.clear();
990    
991      if(!l1)return;
992    
993    
994      // ==========================
995      // CALORIMETER CLUSTERING
996      // ==========================
997      CaloStrip st = CaloStrip(false);
998      ExtHit *caloCl=NULL;
999      Int_t nCl = 0;
1000      for(Int_t ih=0; ih<l1->istrip; ih++){//loop over calorimeter hits
1001        Int_t view = -1;  // 0-x 1-y
1002        Int_t plane = -1; // 0-21
1003        Int_t strip = -1; // 0-95
1004        Float_t mip = l1->DecodeEstrip(ih,view,plane,strip);
1005        
1006        // if(strip == 0 || strip == 95)cout <<" strip "<<strip<<" OK "<<endl;
1007        if(strip < 0 || strip > 95)cout <<" strip "<<strip<<" AHI AHI "<<endl;
1008        
1009        if(view<0)continue;
1010        if(plane<0)continue;
1011        if(strip<0)continue;
1012        if(mip<=0)continue;
1013        
1014        //      cout << ih << " -- "<<view<<" "<<plane<<" "<<strip<<endl;
1015        st.Set(view,plane,strip);// evaluate strip coordinates
1016        //
1017        int   vv = plane*2+(int)(!view);           //0-43 odd x even y
1018        float cc = (view ? st.GetY() : st.GetX()); //coordinate
1019        float ss = mip;                            //signal
1020        float pp = (float)strip;                   //pitch units
1021    
1022        if( vv >= _alg_nViewCal ){//view not included
1023          if(!caloCl)continue;
1024          if(caloCl->mult>0){ //if there is a cluster under construction...
1025            _cal_cl.push_back(ExtHit(*caloCl));//store it
1026            map.insert( pair<int,int>(caloCl->view + vOffset, nCl ));
1027            if(_debug)cout << nCl << " - view(0-43) "<< caloCl->view<<" strip "<<caloCl->start<<"-"<<caloCl->start+caloCl->mult-1<<endl;    
1028            nCl++;      
1029            caloCl->Reset();        
1030            caloCl->Set(strip,vv); // start new one... (without adding the strip)
1031          }
1032          continue;     //... and skip this view
1033        }    
1034        // --------------
1035        // @first hit:
1036        // --------------
1037        if(!caloCl){
1038          caloCl = new ExtHit(strip,vv);//create cluster starting
1039          caloCl->Add(cc,pp,ss);//add strip
1040          continue;
1041        }
1042        //-------------------
1043        //check previous hit
1044        //-------------------    
1045        if(          
1046           vv != caloCl->view ||                  //hit on a different view
1047           strip > caloCl->start + caloCl->mult ||//not adjacent to previous one
1048           strip == 32 || //different ladder
1049           strip == 64 || //different ladder
1050           false ){      
1051          // is a new cluster:
1052          // store previous one, if not already stored
1053          if(caloCl->mult > 0){
1054            _cal_cl.push_back(ExtHit(*caloCl));
1055            map.insert( pair<int,int>(caloCl->view + vOffset, nCl ));
1056            if(_debug)cout << nCl << " - view(0-43) "<< caloCl->view<<" strip "<<caloCl->start<<"-"<<caloCl->start+caloCl->mult-1<<endl;    
1057            nCl++;      
1058          };
1059          //
1060          caloCl->Reset();      // reset
1061          caloCl->Set(strip,vv); // start new one from strip
1062          
1063        }
1064        caloCl->Add(cc,pp,ss); //add strip    
1065      };//end loop over calorimeter hits
1066    
1067      if(caloCl)delete caloCl;
1068    
1069      return;
1070    
1071    };
1072    /**
1073     * Evaluate positions
1074     */
1075    bool ExtTrkingAlg::EvaluateClusterPosition_Tracker( int iclx, int icly,int lad,  int sen, float *xmABar, float *ymABar, float *zmAB){
1076    
1077      TrkParams::Load(1);
1078      if( !TrkParams::IsLoaded(1) ){
1079        cout << "Bool_t ExtTrkingAlg::EvaluateClusterPosition_Tracker ---ERROR--- m.field not loaded "<<endl;
1080        return false;
1081      }    
1082      TrkParams::Load(4);
1083      if( !TrkParams::IsLoaded(4) ){
1084        cout << "Bool_t ExtTrkingAlg::EvaluateClusterPosition_Tracker ---ERROR--- p.f.a. par. not loaded "<<endl;
1085        return false;
1086      }
1087      TrkParams::Load(5);
1088      if( !TrkParams::IsLoaded(5) ){
1089        cout << "Bool_t ExtTrkingAlg::EvaluateClusterPosition_Tracker ---ERROR--- alignment par. not loaded "<<endl;
1090        return false;
1091      }
1092      //    cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;;
1093      int icx = iclx+1;//0=no-cluster,1-N
1094      int icy = icly+1;//0=no-cluster,1-N
1095      int sensor = sen+1;//<< convenzione "Paolo"
1096      int ip = (_trk_l1->GetCluster(TMath::Max(iclx,icly))->view-1)/2;//ip=0-5
1097      if(ip==5 && sensor!=0)sensor=3-sensor;//<< convenzione "Elena"
1098      int ladder = lad+1;
1099      float ax = (float)xmABar[3];
1100      float ay = (float)ymABar[3];
1101      float v[3];
1102      v[0]=xmABar[0];
1103      v[1]=ymABar[0];
1104      v[2]=zmAB[0];
1105      float bfx = 10*TrkParams::GetBX(v);//Tesla
1106      float bfy = 10*TrkParams::GetBY(v);//Tesla
1107      int ipp=ip+1;
1108    
1109      //  cout<<"xyzpam_("<<ipp<<","<<icx<<","<<icy<<","<<ladder<<","<<sensor<<"...)"<<endl;
1110      // e` una cosa bruttissima lo so ma che ci posso fare? mica posso fare tutto io!!
1111      // ma lo sai a che ora me so` svegliata stamattina??
1112      xyzpam_(&ipp,&icx,&icy,&ladder,&sensor,&ax,&ay,&bfx,&bfy);
1113      //  cout<<" ip l s------------------- "<<ip<<ladder<<sensor<<" zm "<<track_.zm[ip]<<endl;
1114      extern cMini2track track_;
1115      
1116      xmABar[0]=track_.xm[ip];
1117      ymABar[0]=track_.ym[ip];
1118      zmAB[0]=track_.zm[ip];
1119      
1120      xmABar[1]=track_.xm_a[ip];
1121      ymABar[1]=track_.ym_a[ip];
1122      zmAB[1]=track_.zm_a[ip];//
1123    
1124      
1125      xmABar[2]=track_.xm_b[ip];
1126      ymABar[2]=track_.ym_b[ip];
1127      zmAB[2]=track_.zm_b[ip];//
1128    
1129      xmABar[4]=track_.resx[ip];
1130      ymABar[4]=track_.resy[ip];
1131    
1132      return true;
1133    
1134    }
1135    // bool ExtTrkingAlg::EvaluateClusterPosition_Calorimeter( int icl,  int sensor, float *xmABar, float *ymABar, float *zmAB, float def){
1136    
1137    //   ExtHit cluster = _cal_cl[icl];
1138    //   //  int plane = 6 + cluster.view;            //0...21
1139    //   int plane = (int)(cluster.view / 2 );    //0...21
1140    //   int ladder = (int)(cluster.coordPU / 32);//0...2
1141    //   bool isX = (cluster.view)%2;
1142    //   bool isY = !isX;
1143      
1144    
1145    //   bool ROTO = (sensor>=0);//if sensor is assigned (0,1,2) the strip can be rototraslated
1146    
1147    
1148    //   CaloStrip st = CaloStrip(false);
1149    //   //----------------------------
1150    //   // Z-coordinate
1151    //   //----------------------------
1152    //   st.Set((int)isY,plane,ladder*32);//
1153    
1154    //   float coordCm_Z0 = st.GetZ();
1155      
1156    //   float size = 80.5; //(cm) larghezza del sensore + zone morte
1157    //   //----------------------------
1158    //   // strip-segment coordinates
1159    //   //----------------------------
1160    //   float coordCm_LA = cluster.coordCm;
1161    //   float coordCm_LB = cluster.coordCm;
1162    //   float coordCm_SA = (-1)*size*1.5; //full calorimeter size
1163    //   float coordCm_SB = (+1)*size*1.5; //full calorimeter size
1164    //   float coordCm_ZA = coordCm_Z0;
1165    //   float coordCm_ZB = coordCm_Z0;
1166    
1167    //   //----------------------------
1168    //   // rototraslation
1169    //   //----------------------------
1170    //   if(ROTO){
1171    //     //-----------------------------
1172    //     // sensor origin coordinates
1173    //     //-----------------------------
1174    //     st.Set(   (int)isY, plane, ladder*32);//
1175    //     float coordCm_L0 = ( isX ? st.GetX() : st.GetY());  
1176    //     st.Set((int)(!isY), plane, sensor*32 );
1177    //     float coordCm_S0 = (!isX ? st.GetX() : st.GetY());
1178        
1179    //     coordCm_SA = (-1)*size*0.5 + (sensor-1)*size; //sensor size
1180    //     coordCm_SB = (+1)*size*0.5 + (sensor-1)*size; //sensor size
1181        
1182    //     //-----------------------------
1183    //     // alignment pamameters
1184    //     //-----------------------------
1185    //     float fPitch = 1.;
1186    //     float shift[]= {0.,0.,0.};//(L,S,Z)
1187    //     float alpha  = 0.;
1188    //     //-----------------------------
1189    
1190    //     //coordinate relative al sensore
1191    //     coordCm_LA = coordCm_LA - coordCm_L0;
1192    //     coordCm_SA = coordCm_SA - coordCm_S0;
1193    //     coordCm_ZA = coordCm_ZA - coordCm_Z0;
1194    //     coordCm_LB = coordCm_LB - coordCm_L0;
1195    //     coordCm_SB = coordCm_SB - coordCm_S0;
1196    //     coordCm_ZB = coordCm_ZB - coordCm_Z0;
1197    
1198    //     //variazione di pitch
1199    //     coordCm_LA = coordCm_LA * fPitch;
1200    //     coordCm_LB = coordCm_LB * fPitch;
1201    
1202    //     //rotazione
1203    //     float LA = coordCm_LA - alpha * coordCm_SA;
1204    //     float SA = coordCm_SA + alpha * coordCm_LA;
1205    //     float LB = coordCm_LB - alpha * coordCm_SB;
1206    //     float SB = coordCm_SB + alpha * coordCm_LB;
1207    
1208    //     //traslazione
1209    //     coordCm_LA = LA         + shift[0];
1210    //     coordCm_SA = SA         + shift[1];
1211    //     coordCm_ZA = coordCm_ZA + shift[2];
1212    //     coordCm_LB = LB         + shift[0];
1213    //     coordCm_SB = SB         + shift[1];
1214    //     coordCm_ZB = coordCm_ZB + shift[2];
1215    
1216    //     //sistema di riferimento generale
1217        
1218    //     coordCm_LA = coordCm_LA + coordCm_L0;
1219    //     coordCm_SA = coordCm_SA + coordCm_S0;
1220    //     coordCm_ZA = coordCm_ZA + coordCm_Z0;
1221    //     coordCm_LB = coordCm_LB + coordCm_L0;
1222    //     coordCm_SB = coordCm_SB + coordCm_S0;
1223    //     coordCm_ZB = coordCm_ZB + coordCm_Z0;
1224    
1225    //   }
1226    
1227    //   xmABar[0] = -100.;
1228    //   ymABar[0] = -100.;
1229    //   zmAB[0]   = 0.5*coordCm_ZA+0.5*coordCm_ZB;
1230      
1231    //   xmABar[1] = (isX ? coordCm_LA : coordCm_SA);
1232    //   ymABar[1] = (isY ? coordCm_LA : coordCm_SA);
1233    //   zmAB[1]   = coordCm_ZA;
1234    
1235      
1236    //   xmABar[2]= (isX ? coordCm_LB : coordCm_SB);
1237    //   ymABar[2]= (isY ? coordCm_LB : coordCm_SB);
1238    //   zmAB[2]  = coordCm_ZB;
1239    
1240    
1241    //   //----------------------------
1242    //   // evaluate spatial resolution
1243    //   //----------------------------
1244    
1245    
1246    //   float res_dig = 0.24/sqrt(12.);
1247    
1248    //   float res_ms = 0.;
1249      
1250    //   if(def!=0.){
1251        
1252    //     float rig = 1/def; if(rig<0)rig = -rig;
1253    //     float factor = 0.;
1254    //     factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*xmABar[3]/180.),2);
1255    //     factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*ymABar[3]/180.),2);
1256    //     factor += 1.;
1257    //     factor = TMath::Sqrt(factor);
1258        
1259    //     int view = cluster.view; //0... 43
1260    
1261    //     int nW = (int)((view + 1)/2);
1262    //     float dW = 0.74*factor;//X0
1263    //     float dW_cm = 0.26*factor;//cm
1264    //     float dSi_cm = 0.; //cm
1265    //     if(isY && plane>0) dSi_cm = ( plane%2 ? 0.228 : 0.428);
1266    //     dSi_cm *= factor;
1267    //     // multiple scattering angle
1268    //     float theta = 0.;//ms
1269    //     theta = 0.0136/rig;
1270    //     theta = theta * sqrt(dW) * (1.+0.038*log(dW));
1271    //     float Vt = theta*theta;            //cov(theta,theta)
1272    //     float Vy = Vt * dW_cm * dW_cm /3.; //cov(y,y)
1273    //     float Vty = 0.87 * Vt * Vy;        //cov(theta,y)
1274    //     float Vc = Vy + dSi_cm * dSi_cm * Vt + 2. * dSi_cm * Vty;
1275    //     float rms = sqrt( Vc );    
1276    //     // if(isY)t_track.resy[ip] = sqrt(rms*rms + t_track.resy[ip]*t_track.resy[ip]);
1277    //     // else   t_track.resx[ip] = sqrt(rms*rms + t_track.resx[ip]*t_track.resx[ip]);
1278    //     // if(isY)          VMS[nW] = rms*rms;
1279    //     // if(_debug)cout <<endl<<view<<" nW "<<nW<<" res(MS) "<<rms<<" factor "<<factor;
1280          
1281    //     if(view>=3)cout << " --> NB error computation not accurate ";
1282    //     res_ms = rms;
1283    
1284    //   }
1285    
1286    
1287    
1288    
1289    //   xmABar[4]= (isX ? sqrt(res_dig*res_dig+res_ms*res_ms) : 1000.);
1290    //   ymABar[4]= (isY ? sqrt(res_dig*res_dig+res_ms*res_ms) : 1000.);
1291    
1292    
1293    
1294    // }
1295    bool ExtTrkingAlg::EvaluateClusterPosition_Calorimeter( int icl, int& sisensor, float *xmABar, float *ymABar, float *zmAB, float def){
1296    
1297      ExtHit cluster = _cal_cl[icl];
1298      //  int plane = 6 + cluster.view;            //0...21
1299      int plane = (int)(cluster.view / 2 );    //0...21
1300      int ladder = (int)(cluster.coordPU / 32);//0...2
1301      //  int view = cluster.view;//0...43    // EM, unused variable
1302      bool isX = (cluster.view)%2;
1303      bool isY = !isX;
1304      
1305    
1306    //  int sisensor = -1;
1307    
1308      
1309    //  CaloStripRoto temp;
1310      
1311      bool ROTO = ( xmABar[0]!=0. && ymABar[0]!=0.);
1312    // ------------------------------------------------------------------------------
1313    //
1314    // check if the cluster belongs to the hit sensor (xmABar[0],ymABar[0] must be provided)
1315    // evaluate the hit sensor number
1316    //
1317    // ------------------------------------------------------------------------------
1318      if(ROTO){
1319    
1320    
1321          for(int is=0;is<3; is++){
1322              int sis = GetCaloSiSensor((int)isY,ladder,is);
1323    //        cout << "view, ladder, is :"<<(int)isY<<ladder<<is<<" ---> sis "<<sis<<endl;
1324              int index = (plane*2+1-(int)isY)*9 + sis;
1325              CaloStripRoto extst = _caloStripRoto[index];
1326    //        cout << "index "<<index<<endl;
1327    //        cout << " -> view     "<<extst.GetView()<<endl;
1328    //        cout << " -> plane    "<<extst.GetPlane()<<endl;
1329    //        cout << " -> sisensor "<<extst.GetSiSensor()<<endl;
1330    //        cout << " -> ladder   "<<extst.GetLadder(<<" "<<GetCaloLadder(extst.GetView(),extst.GetSiSensor())<<endl;      
1331              bool inside = extst.SensorContains(xmABar[0],ymABar[0]);
1332    //        cout << " inside "<<sis<<" ?? "<<inside<<endl;
1333              if( inside ){
1334                  sisensor = sis;
1335                  continue;
1336              }
1337          }
1338    
1339          if(sisensor==-1){
1340    //        cout << "ATTENZIONE!!! sisensor "<<sisensor<<" x,y "<<xmABar[0]<<" "<<ymABar[0]<<endl;
1341          }
1342      }
1343    
1344    // ------------------------------------------------------------------------------
1345    //
1346    // if no sensor is hit, skip rototraslation
1347    //
1348    // ------------------------------------------------------------------------------
1349      if(sisensor==-1)ROTO=false;
1350    // ------------------------------------------------------------------------------
1351    //
1352    // now evaluate strip A-B coordinates
1353    //
1354    // ------------------------------------------------------------------------------
1355      if(ROTO){
1356    
1357          int index = (plane*2+1-(int)isY)*9 + sisensor;
1358          CaloStripRoto extst = _caloStripRoto[index];
1359    //      CaloStrip st = extst.st;
1360          
1361    //       cout << " extst.Set("<<(int)isY<<", "<<plane<<", "<<cluster.coordPU<<", "<<sensor<<") "<<endl;
1362    //       extst.Set((int)isY, plane, cluster.coordPU, sensor);
1363    
1364          if( sisensor != extst.GetSiSensor() )cout << " extst.GetSiSensor() "<<extst.GetSiSensor()<<" != "<<sisensor<<endl;
1365          if( plane != extst.GetPlane() )cout << " extst.GetPlane() "<<extst.GetPlane()<<" != "<<plane<<endl;
1366          if( (int)isY != extst.GetView() )cout << " extst.GetView() "<<extst.GetView()<<" != "<<isY<<endl;
1367    //      cout << " extst.SetStrip("<<cluster.coordPU<<") "<<endl;
1368          extst.SetStrip(cluster.coordPU);
1369          
1370    
1371          xmABar[0] = -100.;
1372          ymABar[0] = -100.;
1373          zmAB[0]   = 0.5*extst.fZA+0.5*extst.fZB;
1374          
1375          xmABar[1] = extst.fXA;
1376          ymABar[1] = extst.fYA;
1377          zmAB[1]   = extst.fZA;
1378          
1379          
1380          xmABar[2]= extst.fXB;
1381          ymABar[2]= extst.fYB;
1382          zmAB[2]  = extst.fZB;
1383    
1384          xmABar[4]= (isX ? extst.GetSpatialResolution(def, xmABar[3], ymABar[3],  1.) : 1000.);
1385          ymABar[4]= (isY ? extst.GetSpatialResolution(def, xmABar[3], ymABar[3],  1.) : 1000.);
1386      
1387    
1388      }else{
1389          
1390    
1391          CaloStrip st = CaloStrip(false);
1392          //----------------------------
1393          // Z-coordinate
1394          //----------------------------
1395          st.Set((int)isY,plane,ladder*32);//
1396    
1397          float coordCm_Z0 = st.GetZ();
1398      
1399          float size = 80.5; //(cm) larghezza del sensore + zone morte
1400          //----------------------------
1401          // strip-segment coordinates
1402          //----------------------------
1403          float coordCm_LA = cluster.coordCm;
1404          float coordCm_LB = cluster.coordCm;
1405          float coordCm_SA = (-1)*size*1.5; //full calorimeter size
1406          float coordCm_SB = (+1)*size*1.5; //full calorimeter size
1407          float coordCm_ZA = coordCm_Z0;
1408          float coordCm_ZB = coordCm_Z0;
1409    
1410          //----------------------------
1411          // rototraslation
1412          //----------------------------
1413    
1414          xmABar[0] = -100.;
1415          ymABar[0] = -100.;
1416          zmAB[0]   = 0.5*coordCm_ZA+0.5*coordCm_ZB;
1417      
1418          xmABar[1] = (isX ? coordCm_LA : coordCm_SA);
1419          ymABar[1] = (isY ? coordCm_LA : coordCm_SA);
1420          zmAB[1]   = coordCm_ZA;
1421    
1422      
1423          xmABar[2]= (isX ? coordCm_LB : coordCm_SB);
1424          ymABar[2]= (isY ? coordCm_LB : coordCm_SB);
1425          zmAB[2]  = coordCm_ZB;
1426    
   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 << "|--|--|--|--|--|--|--|--|--|--|--|"<<endl;  
1427        
1428    
1429    
1430          //----------------------------
1431          // evaluate spatial resolution
1432          //----------------------------
1433    
1434    
1435          float res_dig = 0.24/sqrt(12.);
1436    
1437          float res_ms = 0.;
1438      
1439          if(def!=0.){
1440        
1441              float rig = 1/def; if(rig<0)rig = -rig;
1442              float factor = 0.;
1443              factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*xmABar[3]/180.),2);
1444              factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*ymABar[3]/180.),2);
1445              factor += 1.;
1446              factor = TMath::Sqrt(factor);
1447        
1448              int view = cluster.view; //0... 43
1449    
1450              //      int nW = (int)((view + 1)/2); // EM, unused variable
1451              float dW = 0.74*factor;//X0
1452              float dW_cm = 0.26*factor;//cm
1453              float dSi_cm = 0.; //cm
1454              if(isY && plane>0) dSi_cm = ( plane%2 ? 0.228 : 0.428);
1455              dSi_cm *= factor;
1456              // multiple scattering angle
1457              float theta = 0.;//ms
1458              theta = 0.0136/rig;
1459              theta = theta * sqrt(dW) * (1.+0.038*log(dW));
1460              float Vt = theta*theta;            //cov(theta,theta)
1461              float Vy = Vt * dW_cm * dW_cm /3.; //cov(y,y)
1462              float Vty = 0.87 * Vt * Vy;        //cov(theta,y)
1463              float Vc = Vy + dSi_cm * dSi_cm * Vt + 2. * dSi_cm * Vty;
1464              float rms = sqrt( Vc );        
1465              // if(isY)t_track.resy[ip] = sqrt(rms*rms + t_track.resy[ip]*t_track.resy[ip]);
1466              // else   t_track.resx[ip] = sqrt(rms*rms + t_track.resx[ip]*t_track.resx[ip]);
1467              // if(isY)          VMS[nW] = rms*rms;
1468              // if(_debug)cout <<endl<<view<<" nW "<<nW<<" res(MS) "<<rms<<" factor "<<factor;
1469          
1470              if(view>=3)cout << " --> NB error computation not accurate ";
1471              res_ms = rms;
1472    
1473          }
1474    
1475          xmABar[4]= (isX ? sqrt(res_dig*res_dig+res_ms*res_ms) : 1000.);
1476          ymABar[4]= (isY ? sqrt(res_dig*res_dig+res_ms*res_ms) : 1000.);
1477    
1478      }
1479    
1480      return true; // EM, warning: no return statement in function returning non-void [-Wreturn-type]
1481  }  }
1482    
1483    
1484    
1485    /**
1486     * Process the event
1487     */
1488    void ExtTrkingAlg::ProcessEvent1(Bool_t force){
1489    
1490      if(_debug)cout << " |---------------------------------------------------| "<<endl;
1491      if(_debug)cout << " void ExtTrkingAlg::ProcessEvent1("<<force<<") "<<endl;
1492    
1493      
1494      if(_debug && !_trk_l1)cout << " Missing TrkLevel1 object --- ciao! "<<endl;
1495      if(!_trk_l1)return;
1496    
1497    
1498      //===========================
1499      // pre-selection  
1500      //===========================
1501      
1502      Bool_t RETRACK = CheckEvent();
1503    
1504      if(_debug)cout <<" Pre-selection result: "<<RETRACK<<endl;
1505      if( !force && !RETRACK )return;
1506      
1507    
1508      //===========================
1509      // tracking
1510      //===========================
1511    
1512      // TClonesArray &tracks = *_trkArray;
1513      // tracks.Clear();
1514      TClonesArray &tracks = *_trkArray;
1515      tracks.Clear();
1516    
1517    
1518      multimap<int,ExtTrack> trackCandidates[NEXTVIEWS]; // map: last hit view vs TrkTrack
1519      int mapIndex = 0;
1520      multimap<int,ExtTrack>::iterator trkBegin;
1521      multimap<int,ExtTrack>::iterator trkEnd;
1522      int iCand;
1523    
1524      // ------------------------------------------------
1525      // fill a map VIEW-CLUSTER
1526      // ------------------------------------------------
1527      multimap<int,int> clusterMap;
1528      FillClusterMap(clusterMap,_trk_l1,0);// trk clusters    
1529      FillClusterMap(clusterMap,_cal_l1,12);// calorimeter clusters
1530      
1531    
1532      // ------------------------------------------------
1533      // define iterators over clusters on the same view
1534      // ------------------------------------------------
1535      pair <multimap<int,int>::iterator, multimap<int,int>::iterator> ret[NEXTVIEWS];
1536      multimap<int,int>::iterator it[NEXTVIEWS];
1537      for(int iv=0; iv<NEXTVIEWS; iv++)ret[iv] = clusterMap.equal_range(iv);
1538    
1539      
1540      //====================================================================
1541      // PARAMETERS
1542      //====================================================================
1543      // ------------------------------------------------
1544      // number of points per candidate
1545      // ------------------------------------------------
1546      int nClFixX =  _alg_nClFixX; // n.hits required on X view
1547      int nClFixY = _alg_nClFixY; // n.hits required on X view
1548      int nTrackMAX =  _alg_nTrackMAX; //
1549      //==================================================
1550      // ------------------------------------------------
1551      // start: one candidate for each cluster
1552      // ------------------------------------------------
1553      //==================================================
1554      if(_debug)cout<<"-------- MINIMAL TRACK-CANDIDATES --------";
1555      if(_debug)cout<<nClFixX<<"X"<<nClFixY<<"Y"<<endl;
1556      for(int iv=0; iv<NEXTVIEWS; iv++){                                      //loop over all views
1557        
1558        if( !(iv%2) && (int)(0.5*iv) >  NEXTVIEWS/2-nClFixX)continue;
1559        if(  (iv%2) && (int)(0.5*iv) >  NEXTVIEWS/2-nClFixY)continue;
1560    
1561        if(iv >= NEXTVIEWS - _alg_nViewCal )continue; //skip calo hits
1562        
1563        for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view
1564    
1565    
1566    
1567          int icl = it[iv]->second;
1568    
1569          TrkCluster *cluster = _trk_l1->GetCluster(icl);
1570          //      cluster->Dump();
1571          int ladder = cluster->GetLadder()-1; //
1572          int plane = cluster->GetPlane()-1;
1573          bool isY = (cluster->view)%2;// F77 1-12 ??
1574          bool isX = !isY;
1575          
1576    
1577          /// TRACKER ==============================================================================
1578          for (int is=0; is<2; is++){                                  //loop over sensors
1579    
1580            //cout <<endl<< " view "<<iv<<" plane "<<plane<<" cluster "<<icl<<" on sensor "<<is<<endl;
1581            
1582            ExtTrack& t_track =  *_extTrack;
1583    
1584            t_track.ResetXY();
1585            for(int ip =0 ; ip<t_track.nplanes; ip++ )t_track.SetZ(ip,_zMech[ip]);
1586            //      t_track.Dump();
1587    
1588    
1589            // >>>> insert point
1590            float xmABar[]={0.,0.,0.,0.,0.}; // inizializzare a zero l'angolo, come prima stima per il PFA
1591            float ymABar[]={0.,0.,0.,0.,0.};
1592            float zmAB[]={t_track.zm[plane],0.,0.};//mechanical z-coordinate
1593            EvaluateClusterPosition_Tracker( (isX?icl:-1), (isY?icl:-1), ladder, is, xmABar, ymABar, zmAB);
1594            
1595            if( isX ){
1596              t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
1597              t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
1598              t_track.SetXGood(plane,icl+1,ladder,is);      
1599              //      t_track.SetYGood(plane,0,ladder,is);  
1600            }
1601            if( isY ){
1602              t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
1603              t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
1604              t_track.SetYGood(plane,icl+1,ladder,is);
1605              //      t_track.SetXGood(plane,0,ladder,is);
1606            }
1607            trackCandidates[mapIndex].insert(pair<int,ExtTrack>(iv, t_track));
1608    
1609          }
1610          
1611    
1612    
1613        }
1614      }
1615    
1616      //==============================================================
1617      // -------------------------------------------------------------
1618      // next: add views to each candidate, up to nClFixX+nClFixY
1619      // -------------------------------------------------------------
1620      //==============================================================
1621      for( ;mapIndex < nClFixX+nClFixY-1 ;++mapIndex){
1622    
1623        trkBegin = trackCandidates[mapIndex].begin();
1624        trkEnd   = trackCandidates[mapIndex].end();
1625        if(_debug)cout<<"MAP "<<mapIndex<<" ----------- n.candidates "<<trackCandidates[mapIndex].size()<<endl;
1626    
1627        for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
1628    
1629    
1630          ExtTrack trackCand = (*trkIter).second;
1631          int lastView = (*trkIter).first; //last inserted view
1632          bool lastIsTrk = (lastView < NEXTVIEWS - _alg_nViewCal);
1633          int lastPlane = (int)(0.5*lastView);
1634          if(!lastIsTrk)lastPlane =lastView-6;
1635          int lastLadder = trackCand.GetLadder(lastPlane);
1636          int lastSensor = trackCand.GetSensor(lastPlane);
1637    
1638          
1639          // if(_debug)cout<<"MAP "<<mapIndex<<" last v "<<lastView<<" NX "<<trackCand.GetNX()<<" NY "<<trackCand.GetNY()<<endl;
1640    
1641          for(int iv=lastView+1; iv<NEXTVIEWS; iv++){                             //loop over next views
1642    
1643            bool isTrk = (iv < NEXTVIEWS - _alg_nViewCal);
1644            bool isX = (iv)%2;
1645            bool isY = !isX;
1646            
1647                        
1648            // if(  (iv%2) && (int)(0.5*iv) > 6-nClFixX+trackCand.GetNX())continue; //not enough x-view left
1649            // if( !(iv%2) && (int)(0.5*iv) > 6-nClFixY+trackCand.GetNY())continue; //not enough y-view left
1650    
1651            // if(  (iv%2) && trackCand.GetNX()==nClFixX )continue;//ok?
1652            // if( !(iv%2) && trackCand.GetNY()==nClFixY )continue;//ok?
1653    
1654            if( isX && trackCand.GetNX()==nClFixX )continue;//ok?
1655            if( isY && trackCand.GetNY()==nClFixY )continue;//ok?
1656    
1657    
1658            if( isX && (int)(0.5*iv) >  6+_alg_nViewCal/2     - nClFixX+trackCand.GetNX())continue; //not enough x-view left
1659            if( isY && (int)(0.5*iv) >  6+(_alg_nViewCal+1)/2 - nClFixY+trackCand.GetNY())continue; //not enough y-view left
1660    
1661    
1662            for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view
1663    
1664              int icl = it[iv]->second;
1665    
1666              if(isTrk){
1667                /// TRACKER ==================================================
1668                
1669    
1670                TrkCluster *cluster = _trk_l1->GetCluster(icl);
1671                int ladder = cluster->GetLadder()-1;//(_trk_l1->GetCluster(icl)->maxs-1)%1024 + 1;
1672                int plane = cluster->GetPlane()-1;
1673                // bool isY = (cluster->view)%2;
1674                // bool isX = !isY;
1675                
1676                if( plane==lastPlane && ladder!=lastLadder)continue;
1677                
1678    
1679                for (int is=0; is<2; is++){                                  //loop over sensors
1680                      
1681                  if( plane==lastPlane && is!=lastSensor)continue;
1682                  
1683                  // if(_debug)cout<<"MAP "<<mapIndex<<" last v "<<lastView;
1684                  // if(_debug)cout<<" +cl "<<icl<<" s"<<is<<endl;
1685                  
1686                  ExtTrack t_track = ExtTrack();
1687                  trackCand.Copy(t_track);                             //copy the candidate parameters...
1688                  
1689                  bool isXY = (isX && t_track.YGood(plane)) || (isY && t_track.XGood(plane)) ;
1690                  
1691    
1692                  float xmABar[]={0.,0.,0.,0.,0.}; //importante inizializzare a zero, come prima stima per il PFA
1693                  float ymABar[]={0.,0.,0.,0.,0.};
1694                  float zmAB[]={t_track.zm[plane],0.,0.};
1695    
1696                  if( isX ){
1697                    int iclx = icl;
1698                    int icly = (isXY ? t_track.GetClusterY_ID(plane) : -1 );
1699    
1700                    EvaluateClusterPosition_Tracker( iclx, icly , ladder, is, xmABar, ymABar, zmAB);
1701    
1702                    if(isXY){
1703                      t_track.SetXY(plane,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
1704                      t_track.SetZ(plane,zmAB[0]);
1705                    }else{
1706                      t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
1707                      t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
1708                    }
1709                    t_track.SetXGood(plane,iclx+1,ladder,is);         //...add a point...
1710                    t_track.SetYGood(plane,icly+1,ladder,is);
1711                  }
1712                  //
1713                  if( isY ){
1714                    int iclx = (isXY ? t_track.GetClusterX_ID(plane) : -1 );
1715                    int icly = icl;
1716    
1717                    EvaluateClusterPosition_Tracker( iclx, icly , ladder, is, xmABar, ymABar, zmAB);
1718    
1719                    if(isXY){
1720                      t_track.SetXY(plane,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
1721                      t_track.SetZ(plane,zmAB[0]);
1722                    }else{
1723                      t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
1724                      t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
1725                    };
1726                    t_track.SetYGood(plane,icly+1,ladder,is);
1727                    t_track.SetXGood(plane,iclx+1,ladder,is);
1728                  }
1729    
1730                  trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(iv, t_track)); //...and store a new candidate
1731                };//end loop over sensors
1732                
1733              }else{              
1734                /// CALORIMETER =================================================
1735                ExtHit cluster = _cal_cl[icl];
1736                int plane = 6 + cluster.view;
1737    
1738                if( plane==lastPlane )continue;
1739    
1740                int ladder = (int)(cluster.coordPU / 32);
1741                int sensor = -1; //not yet known
1742                bool isX = (cluster.view)%2;
1743                bool isY = !isX;
1744    
1745                ExtTrack t_track = ExtTrack();
1746                trackCand.Copy(t_track);  
1747    
1748                float xmABar[]={0.,0.,0.,0.,0.}; //importante inizializzare a zero, come prima stima per il PFA
1749                float ymABar[]={0.,0.,0.,0.,0.}; // sensor is not evaluated
1750                float zmAB[]={0.,0.,0.};
1751                EvaluateClusterPosition_Calorimeter( icl, sensor, xmABar, ymABar, zmAB);
1752    
1753    
1754                if( isX ){
1755                  t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
1756                  t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
1757                  t_track.SetXGood(plane,icl+1, ladder, sensor);         //...add a point.
1758                  t_track.SetYGood(plane,    0, ladder, sensor);     //...add a point...
1759                }      
1760                if( isY ){
1761                  t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
1762                  t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
1763                  t_track.SetYGood(plane,icl+1, ladder, sensor);
1764                  t_track.SetXGood(plane,    0, ladder, sensor);
1765                }
1766                trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(iv, t_track)); //...and store a new ca
1767    
1768              }
1769    
1770              if( trackCandidates[mapIndex].size() > (uint)nTrackMAX ){ // EM, compilation warning comparison between signed and unsigned
1771                if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<" > "<<nTrackMAX<<endl;
1772                return;//to many candidates
1773              }
1774    
1775            };//end loop over clusters on the view iv
1776          };//endl loop over views
1777        }//end loop over candidates
1778      }//end iterations
1779      
1780    
1781      // if( trackCandidates[mapIndex].size() > nTrakMAX ){
1782      //   if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<" > "<<nTrackMAX<<endl;
1783      //   return;//to many candidates
1784      // }
1785    
1786      // if(_debug){
1787      //   trkBegin = trackCandidates[mapIndex].begin();
1788      //   trkEnd   = trackCandidates[mapIndex].end();
1789      //   for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
1790      //     ExtTrack &trackCand = (*trkIter).second;
1791      //     cout << " >> ";
1792      //     cout << " X: ";
1793      //     for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
1794      //     cout << " Y: ";
1795      //     for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
1796            
1797      //     cout << endl;
1798      //   }
1799      // }
1800    
1801      if(_debug)cout<<"MAP "<<mapIndex<<" ----------- n.candidates "<<trackCandidates[mapIndex].size()<<endl;
1802    
1803      if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
1804    
1805      if(trackCandidates[mapIndex].size()==0)return;//no good candidates
1806    
1807      //===============================================================================================
1808      // -------------------------------------------------------------
1809      // first selection of track candidates, based on chi2
1810      // -------------------------------------------------------------
1811      //===============================================================================================
1812      if(_debug)cout<<"-------- CHI2 --------"<<endl;
1813      trkBegin = trackCandidates[mapIndex].begin();
1814      trkEnd   = trackCandidates[mapIndex].end();
1815      mapIndex++;
1816      iCand = 0;
1817      for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
1818        int lastView = (*trkIter).first;
1819        ExtTrack &trackCand = (*trkIter).second;      
1820        int ifail=0;
1821        trackCand.ResetFit();
1822        trackCand.Fit(0.,ifail,0);
1823        if(ifail!=0)trackCand.ResetFit();
1824    
1825        // -----------------------------------
1826        // first selection of track candidates
1827        // -----------------------------------
1828        if(TMath::IsNaN(trackCand.chi2))continue;
1829        if(trackCand.chi2 <= 0. || trackCand.chi2 > 1.e7)continue;
1830        //    if(trackCand.chi2 < 0. || trackCand.chi2 > 10. )continue;
1831        iCand++;    
1832              
1833        if(_debug){    
1834         cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
1835          cout <<" | al: ";
1836          for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
1837          cout <<" | ";
1838          cout << " X: ";
1839          for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
1840          // for(int ip=0; ip<trackCand.nplanes; ip++)
1841          // if(trackCand.GetClusterX_ID(ip)>=0)cout << trackCand.GetClusterX_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";
1842          cout << " | ";
1843          cout << " Y: ";
1844          for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
1845          cout << endl;
1846          // for(int ip=0; ip<trackCand.nplanes; ip++)
1847          //        if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";       cout << endl;
1848        }
1849    
1850        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1851        // evaluated coordinates (to define GF)
1852        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1853        // int    ngf = TrkParams::nGF;
1854        // float *zgf = TrkParams::zGF;
1855        // Trajectory tgf = Trajectory(ngf,zgf);
1856        // _tgf->DoTrack(trackCand.al);//<<<< integrate the trajectory
1857        // for(int ip=0; ip<ngf; ip++){
1858        //   trackCand.xGF[ip] = tgf->x[ip];
1859        //   trackCand.yGF[ip] = tgf->y[ip];
1860        // }
1861            
1862        trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, trackCand));
1863    
1864      }
1865      if(trackCandidates[mapIndex].size()==0)return;//no good candidates
1866    
1867      if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
1868      //=================================================================
1869      // ----------------------------------------------------------------
1870      // CALORIMETER
1871      // ----------------------------------------------------------------
1872      //=================================================================
1873      if(_cal_l1){
1874        if(_debug)cout<<"-------- CALORIMETER --------"<<endl;
1875        // -----------------------------------
1876        // evaluate calorimeter variables
1877        // -----------------------------------
1878        // float caloCoord[2][22][96];
1879        // float caloZ[44];
1880        // for(int view=0; view<2; view++){
1881        //   for(int plane=0; plane<22; plane++){
1882        //  for(int strip=0; strip<96; strip++){
1883        //    CaloStrip st = CaloStrip(false);
1884        //    st.Set(view,plane,strip);
1885        //    caloCoord[view][plane][strip]=(view ? st.GetY() : st.GetX());
1886        //    caloZ[plane*2+(view?0:1)]=st.GetZ();
1887        //  }
1888        //   }
1889        // }
1890        // Trajectory caloTj = Trajectory(44,caloZ);
1891        //  for(int iz=0; iz<44;iz++)cout<<caloZ[iz]<<endl;
1892        //test calo
1893        //    float caloEvent[2][22][96];
1894        CaloLevel1* l1=_cal_l1;//event->GetCaloLevel1();
1895        //    if(!l1)return;
1896        vector<float> caloChi2;
1897        vector<int> caloHit; int caloHitMAX=0;
1898        vector<float> caloMip;
1899        trkBegin = trackCandidates[mapIndex].begin();
1900        trkEnd   = trackCandidates[mapIndex].end();
1901        for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){        
1902          ExtTrack trackCand = ((*trkIter).second);  
1903                
1904          // trackCand.DoTrack(&caloTj);//integrate the trajectory through the calorimeter
1905          _caloTj->DoTrack(trackCand.al,trackCand.zini);
1906          //      _caloTj->Dump();
1907          
1908          float chi2 = 0.;
1909          int  nhit = 0;
1910          float cmip = 0.;
1911          for(Int_t ih=0; ih<l1->istrip; ih++){
1912            Int_t view = -1;
1913            Int_t plane = -1;
1914            Int_t strip = -1;
1915            Float_t mip = l1->DecodeEstrip(ih,view,plane,strip);
1916            if(view<0)continue;
1917            if(plane<0)continue;
1918            if(strip<0)continue;
1919            if(mip<=0)continue;
1920            float dxy = _caloCoord[view][plane][strip]-(view?_caloTj->y[plane*2+(view?0:1)]:_caloTj->x[plane*2+(view?0:1)]);
1921            if(plane<11) chi2+= mip * TMath::Power(dxy,2.);
1922            //            caloZ[plane*2+(view?0:1)]=st.GetZ();
1923            if( fabs(dxy) < 2. && plane<11)nhit++;
1924            if( fabs(dxy) < 2. && plane<11)cmip+=mip;
1925          };
1926          if(l1->istrip>0)chi2 = chi2/l1->istrip;
1927                
1928          caloHit.push_back(nhit);
1929          if(nhit>caloHitMAX)caloHitMAX=nhit;
1930    
1931          if(_debug){          
1932            cout << " CALO ";
1933            //      cout << " chi2 "<< chi2 <<" ";
1934            cout << " nhit (< 11th plane) "<< nhit <<" ";
1935            //      cout << " mip  "<< cmip <<" ";
1936            cout <<endl;
1937          };
1938        }
1939    
1940              
1941        //=================================================================
1942        // ----------------------------------------------------------------
1943        // selection of track candidates, based on n.calorimeter hits
1944        // ----------------------------------------------------------------
1945        //=================================================================
1946        trkBegin = trackCandidates[mapIndex].begin();
1947        trkEnd   = trackCandidates[mapIndex].end();
1948        mapIndex++;
1949        iCand = 0;
1950        for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
1951          int lastView = (*trkIter).first;
1952          ExtTrack &trackCand = (*trkIter).second;
1953          if(caloHit[iCand]==caloHitMAX)trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, trackCand));
1954          iCand++;
1955        }
1956        if(trackCandidates[mapIndex].size()==0)return;//no good candidates
1957        if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
1958    
1959      }
1960    
1961    
1962      //=================================================================
1963      // ----------------------------------------------------------------
1964      // TOF
1965      // ----------------------------------------------------------------
1966      //=================================================================
1967      if(_tof_l2){
1968        if(_debug)cout<<"-------- TOF --------"<<endl;
1969        trkBegin = trackCandidates[mapIndex].begin();
1970        trkEnd   = trackCandidates[mapIndex].end();
1971        ToFLevel2* tl2=_tof_l2;//event->GetToFLevel2();
1972        //   if(!tl2)return;
1973        iCand =0;
1974        vector<int> tofHit; int tofHitMAX=0;
1975        for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
1976          //            int lastView = (*trkIter).first;
1977          ExtTrack &trackCand = (*trkIter).second;
1978    
1979          _tgf->DoTrack(trackCand.al);//<<<< integrate the trajectory
1980    
1981           for(int ip=0; ip<TrkParams::nGF; ip++){
1982            trackCand.xGF[ip] = _tgf->x[ip];
1983            trackCand.yGF[ip] = _tgf->y[ip];
1984          }
1985    
1986          int nhit =0;
1987          for(int iToF=0; iToF<6; iToF++){
1988                  
1989            int iGF = iToF;
1990            if(iToF>3)iGF = iToF+8;
1991            int iPaddle = tl2->GetPaddleIdOfTrack( trackCand.xGF[iGF], trackCand.yGF[iGF], iToF);
1992            if( tl2->HitPaddle(iToF,iPaddle) )nhit++;
1993                  
1994          }
1995          iCand++;
1996                
1997          tofHit.push_back(nhit);
1998          if(nhit>tofHitMAX)tofHitMAX=nhit;
1999    
2000          if(_debug){
2001            cout << " TOF nhit "<< nhit <<" ";
2002            cout <<endl;
2003          }
2004                
2005        }
2006        //=================================================================
2007        // ----------------------------------------------------------------
2008        // selection of track candidates, based on n.tof hits
2009        // ----------------------------------------------------------------
2010        //=================================================================
2011        trkBegin = trackCandidates[mapIndex].begin();
2012        trkEnd   = trackCandidates[mapIndex].end();
2013        mapIndex++;
2014        iCand = 0;
2015        for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
2016          int lastView = (*trkIter).first;
2017          ExtTrack &trackCand = (*trkIter).second;
2018          if(tofHit[iCand]==tofHitMAX)trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, ExtTrack(trackCand)));
2019          iCand++;
2020    
2021        }
2022        if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
2023        if(trackCandidates[mapIndex].size()==0)return;//no good candidates
2024    
2025      }
2026    
2027      if(_debug){    
2028        
2029        cout<< " Minimal-track SUMMARY:"<<endl;
2030    
2031        trkBegin = trackCandidates[mapIndex].begin();
2032        trkEnd   = trackCandidates[mapIndex].end();
2033        for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
2034          ExtTrack &trackCand = (*trkIter).second;
2035          cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
2036          cout <<" | al: ";
2037          for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
2038          cout <<" | ";
2039          cout << " X: ";
2040          //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
2041          for(int ip=0; ip<trackCand.nplanes; ip++){
2042            if(trackCand.GetClusterX_ID(ip)>=0)
2043              cout << setw(2)<<trackCand.GetClusterX_ID(ip);
2044            else
2045              cout << "__";
2046          }
2047          cout << " | ";
2048          cout << " Y: ";
2049          //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
2050          for(int ip=0; ip<trackCand.nplanes; ip++){
2051            if(trackCand.GetClusterY_ID(ip)>=0)
2052              cout << setw(2)<<trackCand.GetClusterY_ID(ip);
2053            else
2054              cout << "__";
2055          }
2056          cout << " S:";for(int ip=0; ip<trackCand.nplanes; ip++)cout<< trackCand.GetSensor(ip)+1;
2057          cout << endl;
2058    
2059        }
2060      }
2061      
2062      //=================================================================
2063      // ----------------------------------------------------------------
2064      // TRACK MERGING
2065      // ----------------------------------------------------------------
2066      //=================================================================
2067      if(_debug)cout<<"-------- MERGING --------";
2068      trkBegin = trackCandidates[mapIndex].begin();
2069      trkEnd   = trackCandidates[mapIndex].end();
2070      mapIndex++;
2071      float chi2Min=1e10;
2072      iCand = 0;
2073      for( multimap<int,ExtTrack>::iterator trkIter1 = trkBegin; trkIter1!=trkEnd; ++trkIter1){ //loop over tracks (1)
2074        
2075        ExtTrack &trackCand1 = (*trkIter1).second;// get candidate 1...
2076                
2077        ExtTrack t_track = ExtTrack();// ...create a new ExtTrack instance...
2078        trackCand1.Copy(t_track);     // ...and make a copy
2079                      
2080        bool HBM = false; //has been merged
2081    
2082        for( multimap<int,ExtTrack>::iterator trkIter2 = trkBegin; trkIter2!=trkEnd; ++trkIter2){//loop over tracks (2)
2083                  
2084          ExtTrack &trackCand2 = (*trkIter2).second; // get candidate 2...
2085    
2086          bool CBM = true;
2087          for(int ip=0; ip<NEXTPLANES; ip++){ //loop over planes
2088            if( t_track.XGood(ip)  &&  
2089                trackCand2.XGood(ip) &&  
2090                ( t_track.GetClusterX_ID(ip)!=trackCand2.GetClusterX_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
2091                true )CBM=false; // if both have a x-hit and it is not the same ---> they cannot be merged
2092            if(!CBM)break;      
2093            if( t_track.YGood(ip)  &&
2094                trackCand2.YGood(ip) &&
2095                ( t_track.GetClusterY_ID(ip)!=trackCand2.GetClusterY_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
2096                true )CBM=false; // if both have a y-hit and it is not the same ---> they cannot be merged
2097            if(!CBM)break;        
2098          }
2099          //      if(_debug)cout << " can be merged ? "<<CBM<<endl;
2100          if(!CBM)continue; //go to the next candidate
2101          ///////////
2102          //merging//
2103          ///////////    
2104          // cout <<endl<<" check merging al (1) ";
2105          // for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
2106          // cout <<endl<<"       with    al (2) ";
2107          // for(int i=0; i<5; i++)cout <<setw(10)<< trackCand2.al[i];
2108      
2109          for(int ip=0; ip<NEXTPLANES; ip++){            
2110            if( !(t_track.XGood(ip)) &&  trackCand2.XGood(ip) ){      
2111              int ygood_temp = t_track.ygood[ip];
2112              
2113              if(t_track.YGood(ip)){
2114                t_track.SetXY(ip,trackCand2.xm[ip],t_track.ym[ip],trackCand2.resx[ip],t_track.resy[ip]);
2115                //      cout<<"confronta "<<trackCand2.zm[ip]<<" - "<<t_track.zm[ip]<<" = "<<trackCand2.zm[ip]-t_track.zm[ip]<<endl;
2116              }else{
2117                t_track.SetX(ip,trackCand2.xma[ip],trackCand2.xmb[ip],trackCand2.yma[ip],trackCand2.ymb[ip],trackCand2.resx[ip]);
2118                t_track.SetZ(ip,0.5*(trackCand2.zma[ip]+trackCand2.zmb[ip]));
2119              }
2120              t_track.xgood[ip] = trackCand2.xgood[ip];//assign cluster + sensor
2121              t_track.ygood[ip] = ygood_temp;//assign cluster + sensor
2122    
2123              HBM = true;    
2124            }              
2125            if( !(t_track.YGood(ip)) &&  trackCand2.YGood(ip) ) {    
2126    
2127              int xgood_temp = t_track.xgood[ip];
2128              
2129              if(t_track.XGood(ip)){
2130                t_track.SetXY(ip,t_track.xm[ip],trackCand2.ym[ip],t_track.resx[ip],trackCand2.resy[ip]);    
2131                //      cout<<"confronta "<<trackCand2.zm[ip]<<" - "<<t_track.zm[ip]<<" = "<<trackCand2.zm[ip]-t_track.zm[ip]<<endl;
2132              }else{
2133                t_track.SetY(ip,trackCand2.xma[ip],trackCand2.xmb[ip],trackCand2.yma[ip],trackCand2.ymb[ip],trackCand2.resy[ip]);
2134                t_track.SetZ(ip,0.5*(trackCand2.zma[ip]+trackCand2.zmb[ip]));
2135              }
2136              t_track.ygood[ip] = trackCand2.ygood[ip];//assign cluster + sensor
2137              t_track.xgood[ip] = xgood_temp;//assign cluster + sensor
2138              HBM = true;
2139            }
2140          };        
2141          
2142        }//end loop over tracks (2)
2143    
2144    
2145        // cout << "<<  ";
2146        // for(int ip=0; ip<6; ip++) cout <<setw(4)<< t_track.GetClusterX_ID(ip);
2147        // for(int ip=0; ip<6; ip++) cout <<setw(4)<< t_track.GetClusterY_ID(ip);
2148        // cout << endl;
2149        // cout<< " TRK chi2 "<<setw(10)<<t_track.chi2;
2150        // cout <<" | al ";
2151        // for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
2152        // cout << endl;
2153    
2154        ///////////////////////////////////////////////////////// FIT START
2155        // int ifail=0;
2156        // //           TrkParams::SetDebugMode();
2157        // t_track.FitReset();       //angles set to zero
2158        // t_track.Fit(0.,ifail,0,1);//re-evaluate cluster position from level1
2159        // if( TMath::IsNaN(t_track.chi2) )continue;
2160        // if( ifail != 0 )continue;
2161        // t_track.Fit(0.,ifail,0,1);//evaluate again taking into account the track angle
2162        // if( TMath::IsNaN(t_track.chi2) )continue;
2163        // if( ifail != 0 )continue;
2164        ///////////////////////////////////////////////////////// FIT FINISH
2165        // cout << "Merged: "<<endl;
2166        // cout << "X "; for(int ip=0; ip<6; ip++)cout << t_track.XGood(ip); cout << endl;
2167        // cout << "Y "; for(int ip=0; ip<6; ip++)cout << t_track.YGood(ip); cout << endl;
2168       //
2169        //check if track has been already inserted
2170        //
2171        multimap<int,ExtTrack>::iterator trkB = trackCandidates[mapIndex].begin();
2172        multimap<int,ExtTrack>::iterator trkE = trackCandidates[mapIndex].end();
2173        //
2174        bool ADD = true;
2175        for( multimap<int,ExtTrack>::iterator trkI = trkB; trkI!=trkE; ++trkI){
2176          ExtTrack &trackC = (*trkI).second;
2177          bool SAME=true;
2178          for(int ip=0; ip<NEXTPLANES; ip++){
2179              bool isTrk = (ip < NEXTPLANES - _alg_nViewCal);
2180              if(          t_track.GetClusterX_ID(ip) != trackC.GetClusterX_ID(ip) )SAME=false;
2181              if(          t_track.GetClusterY_ID(ip) != trackC.GetClusterY_ID(ip) )SAME=false;
2182              if( isTrk && t_track.GetSensor(ip)      != trackC.GetSensor(ip)      )SAME=false;                  
2183          }
2184          if(SAME){
2185            ADD=false;
2186            break;
2187          }
2188        }
2189        //    if(_debug)cout <<" ADD ? "<<ADD<<endl;
2190    
2191          //      cout << "-----> INSERT  "<<endl;
2192          if(_debug){
2193    
2194              cout << endl;          
2195              cout<< " TRK chi2 "<<setw(13)<<t_track.chi2;
2196              cout <<" | al: ";
2197              for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
2198              cout <<" | ";
2199              cout << " X: ";
2200              //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.XGood(ip) ? "o" : "-");
2201              for(int ip=0; ip<t_track.nplanes; ip++){
2202                  if(t_track.GetClusterX_ID(ip)>=0)
2203                      cout << setw(2)<<t_track.GetClusterX_ID(ip);
2204                  else
2205                      cout << "__";
2206              }
2207              cout << " | ";
2208              cout << " Y: ";
2209              //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.YGood(ip) ? "o" : "-");
2210              for(int ip=0; ip<t_track.nplanes; ip++){
2211                  if(t_track.GetClusterY_ID(ip)>=0)
2212                      cout << setw(2)<<t_track.GetClusterY_ID(ip);
2213                  else
2214                      cout << "__";
2215              }
2216              cout << " S:";for(int ip=0; ip<t_track.nplanes; ip++)cout<< t_track.GetSensor(ip)+1;
2217              if(HBM)cout << " >> Merged";
2218    
2219          }
2220    
2221    
2222        if(ADD){
2223    
2224    //      cout << " >> re-evaluate...";
2225          // ///////////////////////////////////////////////////////// FIT START
2226          int ifail=0;      
2227          t_track.ResetFit();  
2228          // t_track.Dump();
2229          t_track.Fit(0.,ifail,0);//
2230          if( ifail != 0 || TMath::IsNaN(t_track.chi2) ){
2231            // if(_debug){
2232            //   t_track.ResetFit();      
2233            //   t_track.Fit(0.,ifail,2);//
2234            // }
2235            continue;
2236          }
2237          
2238          // t_track.Fit(0.,ifail,0,1);//evaluate again taking into account the track angle
2239          // if( TMath::IsNaN(t_track.chi2) )continue;
2240          // if( ifail != 0 )continue;
2241          // ///////////////////////////////////////////////////////// FIT FINISH
2242    
2243          
2244    
2245          // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2246          // re-evaluated coordinates
2247          // and
2248          // evaluated other track info
2249          // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
2250          //      float VMS[22];VMS[0]=0.; // EM, 2245:13: warning: variable ?VMS? set but not used [-Wunused-but-set-variable]
2251          for(int ip=0; ip<NEXTPLANES; ip++){
2252    
2253            if(!t_track.XGood(ip)&&!t_track.YGood(ip))continue;
2254            
2255            int iclx = t_track.GetClusterX_ID(ip);
2256            int icly = t_track.GetClusterY_ID(ip);
2257            int ladder = t_track.GetLadder(ip);//ok
2258            int sensor = t_track.GetSensor(ip);// not yet assigned for calo
2259    
2260            bool isTrk = (ip < NEXTPLANES - _alg_nViewCal);  
2261            bool isXY = t_track.XGood(ip)*t_track.YGood(ip);
2262            bool isX  = t_track.XGood(ip) && !isXY;
2263            bool isY  = t_track.YGood(ip) && !isXY;
2264            
2265            //      cout <<endl<< ip << " xgood "<<t_track.xgood[ip]<<" ygood "<<t_track.ygood[ip];
2266            //      cout <<endl<< ip << " iclx "<<iclx<<" icly "<<icly;
2267    
2268            //      continue;
2269            //=========================================
2270            // evaluate angular factor
2271            //=========================================
2272            float factor = 0.;
2273            factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.axv[ip]/180.),2);
2274            factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.ayv[ip]/180.),2);
2275            factor += 1.;
2276            factor = TMath::Sqrt(factor);
2277            //      factor = 1.;///TEMP
2278    
2279    
2280            //=========================================
2281            // re-evaluate  coordinates
2282            //=========================================
2283            float xmABar[]={t_track.xv[ip],0.,0.,t_track.axv[ip],0.};//importante passare gli angoli e la coordinata!!
2284            float ymABar[]={t_track.yv[ip],0.,0.,t_track.ayv[ip],0.};//importante passare gli angoli e la coordinata!!
2285            float zmAB[]={t_track.zv[ip],0.,0.};
2286                    
2287            if(isTrk){
2288                EvaluateClusterPosition_Tracker( iclx, icly, ladder, sensor, xmABar, ymABar, zmAB);          
2289            }else{      
2290                EvaluateClusterPosition_Calorimeter( (int)fmax((float)iclx,(float)icly), sensor, xmABar, ymABar, zmAB, t_track.al[4]);          
2291                // now the hiyt sensor is assigned
2292            }
2293            if(isXY){
2294              t_track.SetXY(ip,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
2295              t_track.SetZ(ip,zmAB[0]);
2296            }
2297            if(isX ){
2298              t_track.SetX( ip,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
2299              t_track.SetZ(ip,0.5*(zmAB[1]+zmAB[2]));
2300            }
2301            if(isY ){
2302              t_track.SetY( ip,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
2303              t_track.SetZ(ip,0.5*(zmAB[1]+zmAB[2]));
2304            }
2305    
2306            t_track.SetXGood(ip,(sensor == -1 ? 0 : iclx+1),ladder,sensor);
2307            t_track.SetYGood(ip,(sensor == -1 ? 0 : icly+1),ladder,sensor);
2308            
2309            //=========================================
2310            // re-evaluate  resolution (calorimeter only)
2311            // taking into account multiple scattering
2312            //=========================================
2313            // if(!isTrk && ip > NEXTPLANES - _alg_nViewCal){
2314            //   int nW = (int)((ip - 6 + 1)/2);
2315            //   float dW = 0.74*factor;//X0
2316            //   float dW_cm = 0.26*factor;//cm
2317            //   float dSi_cm = 0.; //cm
2318            //   if(isY) dSi_cm = ( ip%2 ? 0.228 : 0.428);
2319            //   dSi_cm *= factor;
2320            //   float theta = 0.;
2321            //   if(t_track.GetRigidity()!=0.)theta = 0.0136/t_track.GetRigidity();
2322            //   theta = theta * sqrt(dW) * (1.+0.038*log(dW));
2323            //   float Vt = theta*theta;            //cov(theta,theta)
2324            //   float Vy = Vt * dW_cm * dW_cm /3.; //cov(y,y)
2325            //   float Vty = 0.87 * Vt * Vy;        //cov(theta,y)
2326            //   float Vc = Vy + dSi_cm * dSi_cm * Vt + 2. * dSi_cm * Vty;
2327            //   float rms = sqrt( Vc + (nW>0 ? VMS[nW-1] : 0.) );    
2328            //   if(_debug)cout <<endl<< ip<<" nW "<<nW<<" res(MS) "<<rms;
2329            //   if(isY)t_track.resy[ip] = sqrt(rms*rms + t_track.resy[ip]*t_track.resy[ip]);
2330            //   else   t_track.resx[ip] = sqrt(rms*rms + t_track.resx[ip]*t_track.resx[ip]);
2331            //   if(isY)          VMS[nW] = rms*rms;
2332            // }
2333    
2334            //=========================================
2335            // evaluate other quantities
2336            //=========================================
2337            float smip,multmax;
2338            if(t_track.XGood(ip)){
2339              if(isTrk){
2340                TrkCluster *cluster = _trk_l1->GetCluster(iclx);
2341                float mip = TrkParams::GetMIP(cluster->GetLadder()-1,cluster->view-1);
2342                smip    = cluster->GetSignal()/(mip>0.?mip:1.);
2343                multmax = cluster->GetMultiplicity()*10000+cluster->maxs;
2344              }else{
2345                ExtHit cluster = _cal_cl[iclx];
2346                smip    = cluster.signal;//mip
2347                multmax = cluster.mult*10000+lrint(cluster.coordPU);
2348              }
2349              t_track.dedx_x[ip]   = smip/(factor>0.?factor:1.);
2350              t_track.multmaxx[ip] = multmax;
2351            }
2352            if(t_track.YGood(ip)){
2353              if(isTrk){
2354                TrkCluster *cluster = _trk_l1->GetCluster(icly);
2355                float mip = TrkParams::GetMIP(cluster->GetLadder()-1,cluster->view-1);
2356                smip    = cluster->GetSignal()/(mip>0.?mip:1.);
2357                multmax = cluster->GetMultiplicity()*10000+cluster->maxs;
2358              }else{
2359                ExtHit cluster = _cal_cl[icly];
2360                smip    = cluster.signal;//mip
2361                multmax = cluster.mult*10000+lrint(cluster.coordPU);
2362              }
2363              t_track.dedx_y[ip]   = smip/(factor>0.?factor:1.);
2364              t_track.multmaxy[ip] = multmax;
2365            }
2366          }
2367          ///////////////////////////////////////////////////////// end additional info
2368    
2369    
2370    
2371    
2372          // ///////////////////////////////////////////////////////// RE-FIT START
2373          t_track.Fit(0.,ifail,0);//evaluate again taking into account the track angle
2374          if( TMath::IsNaN(t_track.chi2) )continue;
2375          if( ifail != 0 )continue;
2376          // ///////////////////////////////////////////////////////// RE-FIT FINISH
2377          if(_debug){
2378    //      cout << endl;
2379    //              cout << " >> al: ";for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
2380    //              cout << " chi2: "<<setw(10)<< t_track.chi2;
2381            //      cout << endl;
2382              cout << endl;
2383    //                cout << " >> al: ";for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
2384    //                cout << " chi2: "<<setw(10)<< t_track.chi2;
2385    //                //            cout << endl;
2386              cout<< " TRK chi2 "<<setw(13)<<t_track.chi2;
2387              cout <<" | al: ";
2388              for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
2389              cout <<" | ";
2390              cout << " X: ";
2391              //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.XGood(ip) ? "o" : "-");
2392              for(int ip=0; ip<t_track.nplanes; ip++){
2393                  if(t_track.GetClusterX_ID(ip)>=0)
2394                      cout << setw(2)<<t_track.GetClusterX_ID(ip);
2395                  else
2396                      cout << "__";
2397              }
2398              cout << " | ";
2399              cout << " Y: ";
2400              //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.YGood(ip) ? "o" : "-");
2401              for(int ip=0; ip<t_track.nplanes; ip++){
2402                  if(t_track.GetClusterY_ID(ip)>=0)
2403                      cout << setw(2)<<t_track.GetClusterY_ID(ip);
2404                  else
2405                      cout << "__";
2406              }
2407              cout << " S:";for(int ip=0; ip<t_track.nplanes; ip++)cout<< t_track.GetSensor(ip)+1;
2408              cout << " >> ADD ";
2409          }
2410    
2411         _caloTj->DoTrack(t_track.al,t_track.zini);
2412          for(int ip=0; ip<_caloTj->npoint; ip++){
2413            t_track.xGF[ip] = _caloTj->x[ip];
2414            t_track.yGF[ip] = _caloTj->y[ip];
2415          }
2416    
2417          //              cout << "ADD "<<endl;
2418          trackCandidates[mapIndex].insert(pair<int,ExtTrack>(t_track.GetNX()+t_track.GetNY(), ExtTrack(t_track)));
2419          if(t_track.GetNX()+t_track.GetNY()>5 && t_track.chi2<chi2Min) chi2Min=t_track.chi2;
2420          iCand++;
2421        }//end add-candidate condition
2422                
2423      }//end first loop on candidates
2424      if(_debug)cout <<endl<< "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
2425    
2426    //  cout << "CHI2 min "<<chi2Min<<endl;
2427      //=================================================================
2428      // ----------------------------------------------------------------
2429      // chi2 selection
2430      // ----------------------------------------------------------------
2431      //=================================================================
2432      pair <multimap<int,ExtTrack>::iterator, multimap<int,ExtTrack>::iterator> nxny[NEXTVIEWS];
2433    //  for(int inxny=0; inxny<NEXTVIEWS; inxny++){
2434      bool NMTI =false;//not-minimal track inserted
2435      for(int inxny=NEXTVIEWS-1; inxny>=nClFixX+nClFixY; inxny--){
2436        nxny[inxny] = trackCandidates[mapIndex].equal_range(inxny);
2437    //    cout << " --- n "<<inxny<<endl;
2438        for ( multimap<int,ExtTrack>::iterator it = nxny[inxny].first; it!=nxny[inxny].second; ++it){
2439          ExtTrack &trackC = (*it).second;
2440    //      cout << " NX+NY "<< trackC.GetNX()+trackC.GetNY() <<" chi2 "<<trackC.chi2<< " R(GV) "<<trackC.GetRigidity()<<endl;
2441          //            if( trackC.chi2==0 )continue;
2442          if( (trackC.GetNX()+trackC.GetNY() > 5 && trackC.chi2>chi2Min  ) ||
2443              (trackC.GetNX()+trackC.GetNY() == 5 && NMTI ) ||
2444              false)continue;      
2445    
2446          trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(inxny, ExtTrack(trackC)));
2447    //      cout << " insert "<<endl;
2448    
2449          if( trackC.GetNX()+trackC.GetNY() > 5 )NMTI=true;
2450        }
2451      }
2452      mapIndex++;
2453    
2454      //  if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
2455    
2456      if(trackCandidates[mapIndex].size()==0)return;//no good candidates
2457    
2458      if(_debug){    
2459        
2460        cout<< " Selected-track SUMMARY:"<<endl;
2461    
2462        trkBegin = trackCandidates[mapIndex].begin();
2463        trkEnd   = trackCandidates[mapIndex].end();
2464        for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
2465          ExtTrack &trackCand = (*trkIter).second;
2466          cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
2467          cout <<" | al: ";
2468          for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
2469          cout <<" | ";
2470          cout << " X: ";
2471          //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
2472          for(int ip=0; ip<trackCand.nplanes; ip++){
2473            if(trackCand.GetClusterX_ID(ip)>=0)
2474              cout << setw(2)<<trackCand.GetClusterX_ID(ip);
2475            else
2476              cout << "__";
2477          }
2478          cout << " | ";
2479          cout << " Y: ";
2480          //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
2481          for(int ip=0; ip<trackCand.nplanes; ip++){
2482            if(trackCand.GetClusterY_ID(ip)>=0)
2483              cout << setw(2)<<trackCand.GetClusterY_ID(ip);
2484            else
2485              cout << "__";
2486          }
2487          cout << " S:";for(int ip=0; ip<trackCand.nplanes; ip++)cout<< trackCand.GetSensor(ip)+1;
2488          cout <<"  MDR="<<1./sqrt(trackCand.coval[4][4])<<endl;
2489          cout << endl;
2490          // for(int ip=0; ip<trackCand.nplanes; ip++)
2491          //        if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";       cout << endl;
2492    
2493        }
2494      }
2495      //=================================================================
2496      // ----------------------------------------------------------------
2497      // save tracks
2498      // ----------------------------------------------------------------
2499      //=================================================================
2500      trkBegin = trackCandidates[mapIndex].begin();
2501      trkEnd   = trackCandidates[mapIndex].end();
2502      iCand = 0;
2503      // TClonesArray &tracks = *_trkArray;
2504      // tracks.Clear();
2505      for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){            
2506        ExtTrack trackCand = ((*trkIter).second);
2507        //    trackCand.seqno = iCand;
2508        //    trackCand.Dump();
2509        new(tracks[iCand]) ExtTrack(trackCand);
2510        iCand++;
2511      };
2512    
2513    
2514    
2515    };
2516    /**
2517     *
2518     */
2519    
2520    void ExtTrkingAlg::ProcessEvent2(Bool_t force){
2521    
2522    
2523    //    _debug=true;
2524    
2525      if(_debug)cout << " |---------------------------------------------------| "<<endl;
2526      if(_debug)cout << " void ExtTrkingAlg::ProcessEvent2("<<force<<") "<<endl;
2527    
2528      
2529      if(_debug && !_trk_l1)cout << " Missing TrkLevel1 object --- ciao! "<<endl;
2530      if(!_trk_l1)return;
2531    
2532    
2533      //===========================
2534      // pre-selection  
2535      //===========================
2536      
2537      Bool_t RETRACK = CheckEvent();
2538    
2539      if(_debug)cout <<" Pre-selection result: "<<RETRACK<<endl;
2540      if( !force && !RETRACK )return;
2541    
2542      NEXTVIEWS  = 12 +_alg_nViewCal;
2543      NEXTPLANES = 6  +_alg_nViewCal;
2544    
2545      //===========================
2546      // tracking
2547      //===========================
2548    
2549      // TClonesArray &tracks = *_trkArray;
2550      // tracks.Clear();
2551      TClonesArray &tracks = *_trkArray;
2552      tracks.Clear("C");
2553    
2554      int maxSize = (5+5)*(_alg_nViewCal+1);//boh?!
2555      multimap<int,ExtTrack> trackCandidates[maxSize]; // map: last hit view vs TrkTrack
2556      int mapIndex = 0;
2557      multimap<int,ExtTrack>::iterator trkBegin;
2558      multimap<int,ExtTrack>::iterator trkEnd;
2559      int iCand;
2560    
2561      // ------------------------------------------------
2562      // fill a map VIEW-CLUSTER
2563      // ------------------------------------------------
2564      multimap<int,int> clusterMap;
2565      FillClusterMap(clusterMap,_trk_l1,0);// trk clusters    
2566      FillClusterMap(clusterMap,_cal_l1,12);// calorimeter clusters
2567      
2568    
2569      // ------------------------------------------------
2570      // define iterators over clusters on the same view
2571      // ------------------------------------------------
2572      pair <multimap<int,int>::iterator, multimap<int,int>::iterator> ret[NEXTVIEWS];
2573      multimap<int,int>::iterator it[NEXTVIEWS];
2574      for(int iv=0; iv<NEXTVIEWS; iv++)ret[iv] = clusterMap.equal_range(iv);
2575    
2576    
2577      //====================================================================
2578      // PARAMETERS
2579      //====================================================================
2580      // ------------------------------------------------
2581      // number of points per candidate
2582      // ------------------------------------------------
2583      int nClFixX =  _alg_nClFixX; // n.hits required on X view
2584      int nClFixY = _alg_nClFixY; // n.hits required on X view
2585      int nTrackMAX =  _alg_nTrackMAX; //
2586    
2587    
2588    
2589      /////////////////////////////////////////////////////////////////////////
2590      //
2591      //
2592      // start loop over n.calo planes to be included
2593      //
2594      /////////////////////////////////////////////////////////////////////////
2595    
2596      for(int nViewCal=0; nViewCal<=_alg_nViewCal; nViewCal++){
2597    
2598          NEXTVIEWS  = 12 +nViewCal;
2599          NEXTPLANES = 6  +nViewCal;
2600          
2601          if(_debug)cout<<">>>>>>>>>>> INCLUDE "<<nViewCal<<" CALORIMETER VIEWS "<<endl;
2602    //      cout<<">>>>>>>>>>> INCLUDE "<<nViewCal<<" CALORIMETER VIEWS "<<endl;
2603          //==================================================
2604          // ------------------------------------------------
2605          // start: one candidate for each cluster
2606          // ------------------------------------------------
2607          //==================================================
2608          if(_debug)cout<<"-------- MINIMAL TRACK-CANDIDATES --------";
2609          if(_debug)cout<<nClFixX<<"X"<<nClFixY<<"Y"<<endl;
2610          
2611    
2612          for(int iv=0; iv<NEXTVIEWS; iv++){                                      //loop over all views
2613              
2614              if( !(iv%2) && (int)(0.5*iv) >  NEXTVIEWS/2-nClFixX)continue;
2615              if(  (iv%2) && (int)(0.5*iv) >  NEXTVIEWS/2-nClFixY)continue;
2616              
2617              if(iv >= NEXTVIEWS - nViewCal )continue; //skip calo hits
2618              
2619              for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view
2620                  
2621                  
2622    
2623                  int icl = it[iv]->second;
2624    
2625                  TrkCluster *cluster = _trk_l1->GetCluster(icl);
2626                  //      cluster->Dump();
2627                  int ladder = cluster->GetLadder()-1; //
2628                  int plane = cluster->GetPlane()-1;
2629                  bool isY = (cluster->view)%2;// F77 1-12 ??
2630                  bool isX = !isY;
2631          
2632    
2633                  /// TRACKER ==============================================================================
2634                  for (int is=0; is<2; is++){                                  //loop over sensors
2635    
2636                      //cout <<endl<< " view "<<iv<<" plane "<<plane<<" cluster "<<icl<<" on sensor "<<is<<endl;
2637            
2638                      ExtTrack& t_track =  *_extTrack;
2639    
2640                      t_track.ResetXY();
2641    //                for(int ip =0 ; ip<t_track.nplanes; ip++ )t_track.SetZ(ip,_zMech[ip]);
2642    //                t_track.Dump();
2643    
2644    
2645                      // >>>> insert point
2646                      float xmABar[]={0.,0.,0.,0.,0.}; // inizializzare a zero l'angolo, come prima stima per il PFA
2647                      float ymABar[]={0.,0.,0.,0.,0.};
2648                      float zmAB[]={t_track.zm[plane],0.,0.};//mechanical z-coordinate
2649                      EvaluateClusterPosition_Tracker( (isX?icl:-1), (isY?icl:-1), ladder, is, xmABar, ymABar, zmAB);
2650            
2651                      if( isX ){
2652                          t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
2653                          t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
2654                          t_track.SetXGood(plane,icl+1,ladder,is);  
2655                          //          t_track.SetYGood(plane,0,ladder,is);  
2656                      }
2657                      if( isY ){
2658                          t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
2659                          t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
2660                          t_track.SetYGood(plane,icl+1,ladder,is);
2661                          //          t_track.SetXGood(plane,0,ladder,is);
2662                      }
2663                      trackCandidates[mapIndex].insert(pair<int,ExtTrack>(iv, t_track));
2664    
2665                  }
2666          
2667    
2668    
2669              }
2670          }
2671    
2672    
2673          //==============================================================
2674          // -------------------------------------------------------------
2675          // next: add views to each candidate, up to nClFixX+nClFixY
2676          // -------------------------------------------------------------
2677          //==============================================================
2678    //      for( ;mapIndex < nClFixX+nClFixY-1 ;++mapIndex){
2679          for( int iii=0; iii< nClFixX+nClFixY-1; iii++){
2680    
2681    //        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
2682    
2683              trkBegin = trackCandidates[mapIndex].begin();
2684              trkEnd   = trackCandidates[mapIndex].end();
2685              if(_debug)cout<<"MAP "<<mapIndex<<" ----------- n.candidates "<<trackCandidates[mapIndex].size()<<endl;
2686    
2687              for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
2688    
2689    
2690                  ExtTrack trackCand = (*trkIter).second;
2691                  int lastView = (*trkIter).first; //last inserted view
2692                  bool lastIsTrk = (lastView < NEXTVIEWS - nViewCal);
2693                  int lastPlane = (int)(0.5*lastView);
2694                  if(!lastIsTrk)lastPlane =lastView-6;
2695                  int lastLadder = trackCand.GetLadder(lastPlane);
2696                  int lastSensor = trackCand.GetSensor(lastPlane);
2697    
2698    
2699    //            cout <<" >>> lastView "<<lastView<< " NEXTVIEWS "<<NEXTVIEWS<<endl;
2700    
2701                  for(int iv=lastView+1; iv<NEXTVIEWS; iv++){                             //loop over next views
2702    
2703                      bool isTrk = (iv < NEXTVIEWS - nViewCal);
2704                      bool isX = (iv)%2;
2705                      bool isY = !isX;
2706            
2707                        
2708                      // if(  (iv%2) && (int)(0.5*iv) > 6-nClFixX+trackCand.GetNX())continue; //not enough x-view left
2709                      // if( !(iv%2) && (int)(0.5*iv) > 6-nClFixY+trackCand.GetNY())continue; //not enough y-view left
2710    
2711                      // if(  (iv%2) && trackCand.GetNX()==nClFixX )continue;//ok?
2712                      // if( !(iv%2) && trackCand.GetNY()==nClFixY )continue;//ok?
2713    
2714                      if( isX && trackCand.GetNX()==nClFixX )continue;//ok?
2715                      if( isY && trackCand.GetNY()==nClFixY )continue;//ok?
2716    
2717    
2718                      if( isX && (int)(0.5*iv) >  6+nViewCal/2     - nClFixX+trackCand.GetNX())continue; //not enough x-view left
2719                      if( isY && (int)(0.5*iv) >  6+(nViewCal+1)/2 - nClFixY+trackCand.GetNY())continue; //not enough y-view left
2720    
2721                      for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view
2722    
2723                          int icl = it[iv]->second;
2724    
2725    //                    cout <<" >>> icl "<<icl<<endl;
2726    
2727                          if(isTrk){
2728                              /// TRACKER ==================================================
2729                
2730    
2731                              TrkCluster *cluster = _trk_l1->GetCluster(icl);
2732                              int ladder = cluster->GetLadder()-1;//(_trk_l1->GetCluster(icl)->maxs-1)%1024 + 1;
2733                              int plane = cluster->GetPlane()-1;
2734                              // bool isY = (cluster->view)%2;
2735                              // bool isX = !isY;
2736                
2737                              if( plane==lastPlane && ladder!=lastLadder)continue;
2738                
2739    
2740                              for (int is=0; is<2; is++){                                  //loop over sensors
2741                      
2742                                  if( plane==lastPlane && is!=lastSensor)continue;
2743                  
2744                                  // if(_debug)cout<<"MAP "<<mapIndex<<" last v "<<lastView;
2745                                  // if(_debug)cout<<" +cl "<<icl<<" s"<<is<<endl;
2746                  
2747                                  ExtTrack t_track = ExtTrack();
2748                                  trackCand.Copy(t_track);                             //copy the candidate parameters...
2749                  
2750                                  bool isXY = (isX && t_track.YGood(plane)) || (isY && t_track.XGood(plane)) ;
2751                  
2752    
2753                                  float xmABar[]={0.,0.,0.,0.,0.}; //importante inizializzare a zero, come prima stima per il PFA
2754                                  float ymABar[]={0.,0.,0.,0.,0.};
2755                                  float zmAB[]={t_track.zm[plane],0.,0.};
2756    
2757                                  if( isX ){
2758                                      int iclx = icl;
2759                                      int icly = (isXY ? t_track.GetClusterY_ID(plane) : -1 );
2760    
2761                                      EvaluateClusterPosition_Tracker( iclx, icly , ladder, is, xmABar, ymABar, zmAB);
2762    
2763                                      if(isXY){
2764                                          t_track.SetXY(plane,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
2765                                          t_track.SetZ(plane,zmAB[0]);
2766                                      }else{
2767                                          t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
2768                                          t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
2769                                      }
2770                                      t_track.SetXGood(plane,iclx+1,ladder,is);         //...add a point...
2771                                      t_track.SetYGood(plane,icly+1,ladder,is);
2772                                  }
2773                                  //
2774                                  if( isY ){
2775                                      int iclx = (isXY ? t_track.GetClusterX_ID(plane) : -1 );
2776                                      int icly = icl;
2777    
2778                                      EvaluateClusterPosition_Tracker( iclx, icly , ladder, is, xmABar, ymABar, zmAB);
2779    
2780                                      if(isXY){
2781                                          t_track.SetXY(plane,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
2782                                          t_track.SetZ(plane,zmAB[0]);
2783                                      }else{
2784                                          t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
2785                                          t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
2786                                      };
2787                                      t_track.SetYGood(plane,icly+1,ladder,is);
2788                                      t_track.SetXGood(plane,iclx+1,ladder,is);
2789                                  }
2790    
2791                                  trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(iv, t_track)); //...and store a new candidate
2792                              };//end loop over sensors
2793                
2794                          }else{          
2795                              /// CALORIMETER =================================================
2796                              ExtHit cluster = _cal_cl[icl];
2797                              int plane = 6 + cluster.view;
2798    
2799                              if( plane==lastPlane )continue;
2800    
2801                              int ladder = (int)(cluster.coordPU / 32);
2802                              int sensor = -1; //not yet known
2803                              bool isX = (cluster.view)%2;
2804                              bool isY = !isX;
2805    
2806                              ExtTrack t_track = ExtTrack();
2807                              trackCand.Copy(t_track);  
2808    
2809                              float xmABar[]={0.,0.,0.,0.,0.}; //importante inizializzare a zero, come prima stima per il PFA
2810                              float ymABar[]={0.,0.,0.,0.,0.};
2811                              float zmAB[]={0.,0.,0.};
2812                              EvaluateClusterPosition_Calorimeter( icl, sensor, xmABar, ymABar, zmAB);
2813    
2814    
2815                              if( isX ){
2816                                  t_track.SetX(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
2817                                  t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
2818                                  t_track.SetXGood(plane,icl+1, ladder, sensor);         //...add a point.
2819                                  t_track.SetYGood(plane,    0, ladder, sensor);     //...add a point...
2820                              }        
2821                              if( isY ){
2822                                  t_track.SetY(plane,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
2823                                  t_track.SetZ(plane,0.5*(zmAB[1]+zmAB[2]));
2824                                  t_track.SetYGood(plane,icl+1, ladder, sensor);
2825                                  t_track.SetXGood(plane,    0, ladder, sensor);
2826                              }
2827                              trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(iv, t_track)); //...and store a new ca
2828    
2829                          }
2830    
2831                          if( trackCandidates[mapIndex+1].size() > (uint)nTrackMAX ){ //EM, 2826:50: warning: comparison between signed and unsigned integer expressions [-Wsign-compare]
2832                              if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex+1].size()<<" > "<<nTrackMAX<<endl;
2833                              return;//to many candidates
2834                          }
2835    
2836                      };//end loop over clusters on the view iv
2837                  };//endl loop over views
2838              }//end loop over candidates
2839    
2840              mapIndex++;
2841    
2842          }//end iterations
2843      
2844          if(_debug)cout<<"end of cluster inclusion"<<endl;
2845          if(_debug)cout<<"MAP "<<mapIndex<<" ----------- n.candidates "<<trackCandidates[mapIndex].size()<<endl;
2846    //      if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
2847    
2848    
2849          //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
2850          if(trackCandidates[mapIndex].size()==0){
2851    //        mapIndex++;
2852              continue;//increment one calo plane
2853          }
2854          //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
2855    
2856    
2857          //===============================================================================================
2858          // -------------------------------------------------------------
2859          // first selection of track candidates, based on chi2
2860          // -------------------------------------------------------------
2861          //===============================================================================================
2862          if(_debug)cout<<"-------- CHI2 --------"<<endl;
2863    //      if(mapIndex>0)trackCandidates[mapIndex-1].clear();
2864          trkBegin = trackCandidates[mapIndex].begin();
2865          trkEnd   = trackCandidates[mapIndex].end();
2866          mapIndex++;
2867          iCand = 0;
2868          for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
2869              int lastView = (*trkIter).first;
2870              ExtTrack &trackCand = (*trkIter).second;        
2871              int ifail=0;
2872              
2873    
2874              trackCand.ResetFit();
2875              trackCand.Fit(0.,ifail,0);
2876              if(ifail!=0)trackCand.ResetFit();
2877              
2878    //        trackCand.Dump();
2879    
2880    
2881              // -----------------------------------
2882              // first selection of track candidates
2883              // -----------------------------------
2884              if(TMath::IsNaN(trackCand.chi2))continue;
2885              if(trackCand.chi2 <= 0. || trackCand.chi2 > 1.e7)continue;
2886              //    if(trackCand.chi2 < 0. || trackCand.chi2 > 10. )continue;
2887              iCand++;    
2888              
2889              if(_debug){    
2890                  cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
2891                  cout <<" | al: ";
2892                  for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
2893                  cout <<" | ";
2894                  cout << " X: ";
2895                  for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
2896                  // for(int ip=0; ip<trackCand.nplanes; ip++)
2897                  // if(trackCand.GetClusterX_ID(ip)>=0)cout << trackCand.GetClusterX_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";
2898                  cout << " | ";
2899                  cout << " Y: ";
2900                  for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
2901                  cout << endl;
2902                  // for(int ip=0; ip<trackCand.nplanes; ip++)
2903                  //        if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";       cout << endl;
2904              }
2905              trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, trackCand));
2906    
2907          }
2908    
2909    
2910          //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
2911          if(trackCandidates[mapIndex].size()==0){
2912              mapIndex++;
2913              continue;//increment one calo plane
2914          }
2915          //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
2916    
2917          if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
2918    
2919    
2920    
2921          //=================================================================
2922          // ----------------------------------------------------------------
2923          // CALORIMETER
2924          // ----------------------------------------------------------------
2925          //=================================================================
2926          if(_cal_l1){
2927              if(_debug)cout<<"-------- CALORIMETER --------"<<endl;
2928    
2929              CaloLevel1* l1=_cal_l1;//event->GetCaloLevel1();
2930              //    if(!l1)return;
2931              vector<float> caloChi2;
2932              vector<int> caloHit; int caloHitMAX=0;
2933              vector<float> caloMip;
2934    //        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
2935              trkBegin = trackCandidates[mapIndex].begin();
2936              trkEnd   = trackCandidates[mapIndex].end();
2937              for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){          
2938                  ExtTrack trackCand = ((*trkIter).second);  
2939                
2940                  // trackCand.DoTrack(&caloTj);//integrate the trajectory through the calorimeter
2941                  _caloTj->DoTrack(trackCand.al,trackCand.zini);
2942                  //      _caloTj->Dump();
2943          
2944                  float chi2 = 0.;
2945                  int  nhit = 0;
2946                  float cmip = 0.;
2947                  for(Int_t ih=0; ih<l1->istrip; ih++){
2948                      Int_t view = -1;
2949                      Int_t plane = -1;
2950                      Int_t strip = -1;
2951                      Float_t mip = l1->DecodeEstrip(ih,view,plane,strip);
2952                      if(view<0)continue;
2953                      if(plane<0)continue;
2954                      if(strip<0)continue;
2955                      if(mip<=0)continue;
2956                      float dxy = _caloCoord[view][plane][strip]-(view?_caloTj->y[plane*2+(view?0:1)]:_caloTj->x[plane*2+(view?0:1)]);
2957                      if(plane<11) chi2+= mip * TMath::Power(dxy,2.);
2958                      //          caloZ[plane*2+(view?0:1)]=st.GetZ();
2959                      if( fabs(dxy) < 2. && plane<11)nhit++;
2960                      if( fabs(dxy) < 2. && plane<11)cmip+=mip;
2961                  };
2962                  if(l1->istrip>0)chi2 = chi2/l1->istrip;
2963                
2964                  caloHit.push_back(nhit);
2965                  if(nhit>caloHitMAX)caloHitMAX=nhit;
2966    
2967                  if(_debug){          
2968                      cout << " CALO ";
2969                      //    cout << " chi2 "<< chi2 <<" ";
2970                      cout << " nhit (< 11th plane) "<< nhit <<" ";
2971                      //    cout << " mip  "<< cmip <<" ";
2972                      cout <<endl;
2973                  };
2974              }
2975    
2976              //=================================================================
2977              // ----------------------------------------------------------------
2978              // selection of track candidates, based on n.calorimeter hits
2979              // ----------------------------------------------------------------
2980              //=================================================================
2981    //        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
2982              trkBegin = trackCandidates[mapIndex].begin();
2983              trkEnd   = trackCandidates[mapIndex].end();
2984              mapIndex++;
2985              iCand = 0;
2986              for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
2987                  int lastView = (*trkIter).first;
2988                  ExtTrack &trackCand = (*trkIter).second;
2989                  if(caloHit[iCand]==caloHitMAX)trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, trackCand));
2990                  iCand++;
2991              }
2992    
2993              //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
2994              if(trackCandidates[mapIndex].size()==0){
2995                  mapIndex++;
2996                  continue;//increment one calo plane
2997              }
2998              //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
2999    
3000              if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
3001    
3002          }// end if calo event
3003    
3004    
3005    
3006    
3007          //=================================================================
3008          // ----------------------------------------------------------------
3009          // TOF
3010          // ----------------------------------------------------------------
3011          //=================================================================
3012          if(_tof_l2){
3013              if(_debug)cout<<"-------- TOF --------"<<endl;
3014    //        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
3015              trkBegin = trackCandidates[mapIndex].begin();
3016              trkEnd   = trackCandidates[mapIndex].end();
3017              ToFLevel2* tl2=_tof_l2;//event->GetToFLevel2();
3018              //   if(!tl2)return;
3019              iCand =0;
3020              vector<int> tofHit; int tofHitMAX=0;
3021              for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
3022                  //            int lastView = (*trkIter).first;
3023                  ExtTrack &trackCand = (*trkIter).second;
3024    
3025                  _tgf->DoTrack(trackCand.al);//<<<< integrate the trajectory
3026    
3027                  for(int ip=0; ip<TrkParams::nGF; ip++){
3028                      trackCand.xGF[ip] = _tgf->x[ip];
3029                      trackCand.yGF[ip] = _tgf->y[ip];
3030                  }
3031    
3032                  int nhit =0;
3033                  for(int iToF=0; iToF<6; iToF++){
3034                  
3035                      int iGF = iToF;
3036                      if(iToF>3)iGF = iToF+8;
3037                      int iPaddle = tl2->GetPaddleIdOfTrack( trackCand.xGF[iGF], trackCand.yGF[iGF], iToF);
3038                      if( tl2->HitPaddle(iToF,iPaddle) )nhit++;
3039                  
3040                  }
3041                  iCand++;
3042                
3043                  tofHit.push_back(nhit);
3044                  if(nhit>tofHitMAX)tofHitMAX=nhit;
3045    
3046                  if(_debug){
3047                      cout << " TOF nhit "<< nhit <<" ";
3048                      cout <<endl;
3049                  }
3050                
3051        
3052              }
3053              //=================================================================
3054              // ----------------------------------------------------------------
3055              // selection of track candidates, based on n.tof hits
3056              // ----------------------------------------------------------------
3057              //=================================================================
3058    //        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
3059              trkBegin = trackCandidates[mapIndex].begin();
3060              trkEnd   = trackCandidates[mapIndex].end();
3061              mapIndex++;
3062              iCand = 0;
3063              for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
3064                  int lastView = (*trkIter).first;
3065                  ExtTrack &trackCand = (*trkIter).second;
3066                  if(tofHit[iCand]==tofHitMAX)trackCandidates[mapIndex].insert(pair<int,ExtTrack>(lastView, ExtTrack(trackCand)));
3067                  iCand++;
3068    
3069              }
3070              if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
3071    
3072              //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
3073              if(trackCandidates[mapIndex].size()==0){
3074                  mapIndex++;
3075                  continue;//increment one calo plane
3076              }
3077              //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
3078    
3079          }
3080    
3081          if(_debug){    
3082        
3083              cout<< " Minimal-track SUMMARY:"<<endl;
3084    
3085    //        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
3086              trkBegin = trackCandidates[mapIndex].begin();
3087              trkEnd   = trackCandidates[mapIndex].end();
3088              for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
3089                  ExtTrack &trackCand = (*trkIter).second;
3090                  cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
3091                  cout <<" | al: ";
3092                  for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
3093                  cout <<" | ";
3094                  cout << " X: ";
3095                  //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
3096                  for(int ip=0; ip<trackCand.nplanes; ip++){
3097                      if(trackCand.GetClusterX_ID(ip)>=0)
3098                          cout << setw(2)<<trackCand.GetClusterX_ID(ip);
3099                      else
3100                          cout << "__";
3101                  }
3102                  cout << " | ";
3103                  cout << " Y: ";
3104                  //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
3105                  for(int ip=0; ip<trackCand.nplanes; ip++){
3106                      if(trackCand.GetClusterY_ID(ip)>=0)
3107                          cout << setw(2)<<trackCand.GetClusterY_ID(ip);
3108                      else
3109                          cout << "__";
3110                  }
3111                  cout << " S:";for(int ip=0; ip<trackCand.nplanes; ip++)cout<< trackCand.GetSensor(ip)+1;
3112                  cout << endl;
3113    
3114              }
3115          }
3116          //=================================================================
3117          // ----------------------------------------------------------------
3118          // TRACK MERGING
3119          // ----------------------------------------------------------------
3120          //=================================================================
3121          if(_debug)cout<<"-------- MERGING --------";
3122    //        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
3123          trkBegin = trackCandidates[mapIndex].begin();
3124          trkEnd   = trackCandidates[mapIndex].end();
3125          mapIndex++;
3126          float chi2Min=1e10;
3127          iCand = 0;
3128          for( multimap<int,ExtTrack>::iterator trkIter1 = trkBegin; trkIter1!=trkEnd; ++trkIter1){ //loop over tracks (1)
3129        
3130              ExtTrack &trackCand1 = (*trkIter1).second;// get candidate 1...
3131                
3132              ExtTrack t_track = ExtTrack();// ...create a new ExtTrack instance...
3133              trackCand1.Copy(t_track);     // ...and make a copy
3134                      
3135              bool HBM = false; //has been merged
3136    
3137              for( multimap<int,ExtTrack>::iterator trkIter2 = trkBegin; trkIter2!=trkEnd; ++trkIter2){//loop over tracks (2)
3138                  
3139                  ExtTrack &trackCand2 = (*trkIter2).second; // get candidate 2...
3140    
3141                  bool CBM = true;
3142                  for(int ip=0; ip<NEXTPLANES; ip++){ //loop over planes
3143                      if( t_track.XGood(ip)  &&  
3144                          trackCand2.XGood(ip) &&  
3145                          ( t_track.GetClusterX_ID(ip)!=trackCand2.GetClusterX_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
3146                          true )CBM=false; // if both have a x-hit and it is not the same ---> they cannot be merged
3147                      if(!CBM)break;      
3148                      if( t_track.YGood(ip)  &&
3149                          trackCand2.YGood(ip) &&
3150                          ( t_track.GetClusterY_ID(ip)!=trackCand2.GetClusterY_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
3151                          true )CBM=false; // if both have a y-hit and it is not the same ---> they cannot be merged
3152                      if(!CBM)break;        
3153                  }
3154                  //      if(_debug)cout << " can be merged ? "<<CBM<<endl;
3155                  if(!CBM)continue; //go to the next candidate
3156                  ///////////
3157                  //merging//
3158                  ///////////    
3159                  // cout <<endl<<" check merging al (1) ";
3160                  // for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
3161                  // cout <<endl<<"       with    al (2) ";
3162                  // for(int i=0; i<5; i++)cout <<setw(10)<< trackCand2.al[i];
3163                  for(int ip=0; ip<NEXTPLANES; ip++){            
3164                      if( !(t_track.XGood(ip)) &&  trackCand2.XGood(ip) ){            
3165                          int ygood_temp = t_track.ygood[ip];
3166              
3167                          if(t_track.YGood(ip)){
3168                              t_track.SetXY(ip,trackCand2.xm[ip],t_track.ym[ip],trackCand2.resx[ip],t_track.resy[ip]);
3169                              //        cout<<"confronta "<<trackCand2.zm[ip]<<" - "<<t_track.zm[ip]<<" = "<<trackCand2.zm[ip]-t_track.zm[ip]<<endl;
3170                          }else{
3171                              t_track.SetX(ip,trackCand2.xma[ip],trackCand2.xmb[ip],trackCand2.yma[ip],trackCand2.ymb[ip],trackCand2.resx[ip]);
3172                              t_track.SetZ(ip,0.5*(trackCand2.zma[ip]+trackCand2.zmb[ip]));
3173                          }
3174                          t_track.xgood[ip] = trackCand2.xgood[ip];//assign cluster + sensor
3175                          t_track.ygood[ip] = ygood_temp;//assign cluster + sensor
3176    
3177                          HBM = true;        
3178                      }            
3179                      if( !(t_track.YGood(ip)) &&  trackCand2.YGood(ip) ) {          
3180    
3181                          int xgood_temp = t_track.xgood[ip];
3182              
3183                          if(t_track.XGood(ip)){
3184                              t_track.SetXY(ip,t_track.xm[ip],trackCand2.ym[ip],t_track.resx[ip],trackCand2.resy[ip]);      
3185                              //        cout<<"confronta "<<trackCand2.zm[ip]<<" - "<<t_track.zm[ip]<<" = "<<trackCand2.zm[ip]-t_track.zm[ip]<<endl;
3186                          }else{
3187                              t_track.SetY(ip,trackCand2.xma[ip],trackCand2.xmb[ip],trackCand2.yma[ip],trackCand2.ymb[ip],trackCand2.resy[ip]);
3188                              t_track.SetZ(ip,0.5*(trackCand2.zma[ip]+trackCand2.zmb[ip]));
3189                          }
3190                          t_track.ygood[ip] = trackCand2.ygood[ip];//assign cluster + sensor
3191                          t_track.xgood[ip] = xgood_temp;//assign cluster + sensor
3192                          HBM = true;
3193                      }
3194                  };        
3195          
3196              }//end loop over tracks (2)
3197              //
3198              //check if track has been already inserted
3199              //
3200    //        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
3201              multimap<int,ExtTrack>::iterator trkB = trackCandidates[mapIndex].begin();
3202              multimap<int,ExtTrack>::iterator trkE = trackCandidates[mapIndex].end();
3203              //
3204              bool ADD = true;
3205              for( multimap<int,ExtTrack>::iterator trkI = trkB; trkI!=trkE; ++trkI){
3206                  ExtTrack &trackC = (*trkI).second;
3207                  bool SAME=true;
3208                  for(int ip=0; ip<NEXTPLANES; ip++){
3209                      bool isTrk = (ip < NEXTPLANES - nViewCal);
3210                      if(          t_track.GetClusterX_ID(ip) != trackC.GetClusterX_ID(ip) )SAME=false;
3211                      if(          t_track.GetClusterY_ID(ip) != trackC.GetClusterY_ID(ip) )SAME=false;
3212                      if( isTrk && t_track.GetSensor(ip)      != trackC.GetSensor(ip)      )SAME=false;                  
3213                  }
3214                  if(SAME){
3215                      ADD=false;
3216                      break;
3217                  }
3218              }
3219              //    if(_debug)cout <<" ADD ? "<<ADD<<endl;
3220    
3221                  //      cout << "-----> INSERT  "<<endl;
3222              if(_debug){
3223    
3224                  cout << endl;              
3225    //            cout<< " TRK chi2 "<<setw(13)<<t_track.chi2;
3226                  cout<< " TRK chi2 "<<".............";
3227                  cout <<" | al: ";
3228    //            for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
3229                  for(int i=0; i<5; i++)cout <<"..........";
3230                  cout <<" | ";
3231                  cout << " X: ";
3232                  //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.XGood(ip) ? "o" : "-");
3233                  for(int ip=0; ip<t_track.nplanes; ip++){
3234                      if(t_track.GetClusterX_ID(ip)>=0)
3235                          cout << setw(2)<<t_track.GetClusterX_ID(ip);
3236                      else
3237                          cout << "__";
3238                  }
3239                  cout << " | ";
3240                  cout << " Y: ";
3241                  //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.YGood(ip) ? "o" : "-");
3242                  for(int ip=0; ip<t_track.nplanes; ip++){
3243                      if(t_track.GetClusterY_ID(ip)>=0)
3244                          cout << setw(2)<<t_track.GetClusterY_ID(ip);
3245                      else
3246                          cout << "__";
3247                  }
3248                  cout << " S:";for(int ip=0; ip<t_track.nplanes; ip++)cout<< t_track.GetSensor(ip)+1;
3249                  if(HBM)cout << " >> Merged";
3250    //            cout << endl;
3251              }      
3252    
3253              if(ADD){
3254    //            cout << " >> re-evaluate...";
3255    //                cout << " >> Merged? "<<HBM;
3256    //                cout << "| X:";
3257    //                for(int ip=0; ip<NEXTPLANES; ip++){
3258    //                    if(t_track.GetClusterX_ID(ip)>=0)
3259    //                        cout << setw(2)<<t_track.GetClusterX_ID(ip);
3260    //                    else
3261    //                        cout << "__";
3262    //                }
3263    //                //    for(int ip=0; ip<NEXTPLANES; ip++)cout<< t_track.XGood(ip);
3264    //                cout << "| Y:";
3265    //                for(int ip=0; ip<NEXTPLANES; ip++){
3266    //                    if(t_track.GetClusterY_ID(ip)>=0)
3267    //                        cout << setw(2)<<t_track.GetClusterY_ID(ip);
3268    //                    else
3269    //                        cout << "__";
3270    //                }
3271    //                //    for(int ip=0; ip<NEXTPLANES; ip++)cout<< t_track.YGood(ip);    
3272    //                cout << " S:";for(int ip=0; ip<NEXTPLANES; ip++)cout<< t_track.GetSensor(ip)+1;
3273    //                cout << endl;
3274                  
3275                  // ///////////////////////////////////////////////////////// FIT START
3276                  int ifail=0;      
3277                  t_track.ResetFit();  
3278                  // t_track.Dump();
3279                  t_track.Fit(0.,ifail,0);//
3280                  if( ifail != 0 || TMath::IsNaN(t_track.chi2) ){
3281                      // if(_debug){
3282                      //   t_track.ResetFit();      
3283                      //   t_track.Fit(0.,ifail,2);//
3284                      // }
3285                      continue;
3286                  }
3287          
3288                  // t_track.Fit(0.,ifail,0,1);//evaluate again taking into account the track angle
3289                  // if( TMath::IsNaN(t_track.chi2) )continue;
3290                  // if( ifail != 0 )continue;
3291                  // ///////////////////////////////////////////////////////// FIT FINISH
3292    
3293          
3294    
3295                  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3296                  // re-evaluated coordinates
3297                  // and
3298                  // evaluated other track info
3299                  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
3300                  //float VMS[22];VMS[0]=0.;  //EM, 3295:14: warning: variable ?VMS? set but not used [-Wunused-but-set-variable]
3301                  for(int ip=0; ip<NEXTPLANES; ip++){
3302    
3303                      if(!t_track.XGood(ip)&&!t_track.YGood(ip))continue;
3304            
3305                      int iclx = t_track.GetClusterX_ID(ip);
3306                      int icly = t_track.GetClusterY_ID(ip);
3307                      int ladder = t_track.GetLadder(ip);//ok
3308                      int sensor = t_track.GetSensor(ip);//not assigned yet for calo
3309    
3310                      bool isTrk = (ip < NEXTPLANES - nViewCal);      
3311                      bool isXY = t_track.XGood(ip)*t_track.YGood(ip);
3312                      bool isX  = t_track.XGood(ip) && !isXY;
3313                      bool isY  = t_track.YGood(ip) && !isXY;
3314            
3315                      //    cout <<endl<< ip << " xgood "<<t_track.xgood[ip]<<" ygood "<<t_track.ygood[ip];
3316                      //    cout <<endl<< ip << " iclx "<<iclx<<" icly "<<icly;
3317    
3318                      //    continue;
3319                      //=========================================
3320                      // evaluate angular factor
3321                      //=========================================
3322                      float factor = 0.;
3323                      factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.axv[ip]/180.),2);
3324                      factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.ayv[ip]/180.),2);
3325                      factor += 1.;
3326                      factor = TMath::Sqrt(factor);
3327                      //    factor = 1.;///TEMP
3328    
3329                      //=========================================
3330                      // re-evaluate  coordinates
3331                      //=========================================
3332                      float xmABar[]={t_track.xv[ip],0.,0.,t_track.axv[ip],0.};//importante passare gli angoli!!
3333                      float ymABar[]={t_track.yv[ip],0.,0.,t_track.ayv[ip],0.};//importante passare gli angoli!!
3334                      float zmAB[]={t_track.zv[ip],0.,0.};
3335                    
3336                      if(isTrk){
3337                          EvaluateClusterPosition_Tracker( iclx, icly, ladder, sensor, xmABar, ymABar, zmAB);        
3338                      }else{
3339    //                    cout <<" ------------ EvaluateClusterPosition_Calorimeter ----------------"<<endl;
3340    //                    cout << "IN: "<<(int)fmax((float)iclx,(float)icly)<<endl;
3341    //                    cout << "IN: "<<sensor<<endl;
3342    //                    cout << "IN: "<<xmABar[0]<<endl;
3343    //                    cout << "IN: "<<ymABar[0]<<endl;
3344    //                    cout << "IN: "<<zmAB[0]<<endl;
3345    //                    cout << "IN: "<<t_track.al[4]<<endl;
3346                          EvaluateClusterPosition_Calorimeter( (int)fmax((float)iclx,(float)icly), sensor, xmABar, ymABar, zmAB, t_track.al[4]);
3347    //                    cout << ip << "  SENSOR "<<sensor<<endl;
3348                      }
3349                      if(isXY){
3350                          t_track.SetXY(ip,xmABar[0],ymABar[0],xmABar[4],ymABar[4]);
3351                          t_track.SetZ(ip,zmAB[0]);
3352                      }
3353                      if(isX ){
3354                          t_track.SetX( ip,xmABar[1],xmABar[2],ymABar[1],ymABar[2],xmABar[4]);
3355                          t_track.SetZ(ip,0.5*(zmAB[1]+zmAB[2]));
3356                      }
3357                      if(isY ){
3358                          t_track.SetY( ip,xmABar[1],xmABar[2],ymABar[1],ymABar[2],ymABar[4]);
3359                          t_track.SetZ(ip,0.5*(zmAB[1]+zmAB[2]));
3360                      }
3361    
3362                      t_track.SetXGood(ip,(sensor == -1 ? 0 : iclx+1),ladder,sensor);
3363                      t_track.SetYGood(ip,(sensor == -1 ? 0 : icly+1),ladder,sensor);
3364    //                t_track.SetXGood(ip,iclx+1,ladder,sensor);
3365    //                t_track.SetYGood(ip,icly+1,ladder,sensor);
3366            
3367                      //=========================================
3368                      // evaluate other quantities
3369                      //=========================================
3370                      float smip,multmax;
3371                      if(t_track.XGood(ip)){
3372                          if(isTrk){
3373                              TrkCluster *cluster = _trk_l1->GetCluster(iclx);
3374                              float mip = TrkParams::GetMIP(cluster->GetLadder()-1,cluster->view-1);
3375                              smip    = cluster->GetSignal()/(mip>0.?mip:1.);
3376                              multmax = cluster->GetMultiplicity()*10000+cluster->maxs;
3377                          }else{
3378                              ExtHit cluster = _cal_cl[iclx];
3379                              smip    = cluster.signal;//mip
3380                              multmax = cluster.mult*10000+lrint(cluster.coordPU);
3381                          }
3382                          t_track.dedx_x[ip]   = smip/(factor>0.?factor:1.);
3383                          t_track.multmaxx[ip] = multmax;
3384                      }
3385                      if(t_track.YGood(ip)){
3386                          if(isTrk){
3387                              TrkCluster *cluster = _trk_l1->GetCluster(icly);
3388                              float mip = TrkParams::GetMIP(cluster->GetLadder()-1,cluster->view-1);
3389                              smip    = cluster->GetSignal()/(mip>0.?mip:1.);
3390                              multmax = cluster->GetMultiplicity()*10000+cluster->maxs;
3391                          }else{
3392                              ExtHit cluster = _cal_cl[icly];
3393                              smip    = cluster.signal;//mip
3394                              multmax = cluster.mult*10000+lrint(cluster.coordPU);
3395                          }
3396                          t_track.dedx_y[ip]   = smip/(factor>0.?factor:1.);
3397                          t_track.multmaxy[ip] = multmax;
3398                      }
3399                  }
3400                  ///////////////////////////////////////////////////////// end additional info
3401    
3402                  // ///////////////////////////////////////////////////////// RE-FIT START
3403                  t_track.Fit(0.,ifail,0);//evaluate again taking into account the track angle
3404                  if( TMath::IsNaN(t_track.chi2) )continue;
3405                  if( ifail != 0 )continue;
3406                  // ///////////////////////////////////////////////////////// RE-FIT FINISH
3407                  if(_debug){
3408                      cout << endl;
3409    //                cout << " >> al: ";for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
3410    //                cout << " chi2: "<<setw(10)<< t_track.chi2;
3411    //                //            cout << endl;
3412                      cout<< " TRK chi2 "<<setw(13)<<t_track.chi2;
3413                      cout <<" | al: ";
3414                      for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
3415                      cout <<" | ";
3416                      cout << " X: ";
3417                      //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.XGood(ip) ? "o" : "-");
3418                      for(int ip=0; ip<t_track.nplanes; ip++){
3419                          if(t_track.GetClusterX_ID(ip)>=0)
3420                              cout << setw(2)<<t_track.GetClusterX_ID(ip);
3421                          else
3422                              cout << "__";
3423                      }
3424                      cout << " | ";
3425                      cout << " Y: ";
3426                      //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (t_track.YGood(ip) ? "o" : "-");
3427                      for(int ip=0; ip<t_track.nplanes; ip++){
3428                          if(t_track.GetClusterY_ID(ip)>=0)
3429                              cout << setw(2)<<t_track.GetClusterY_ID(ip);
3430                          else
3431                              cout << "__";
3432                      }
3433                      cout << " S:";for(int ip=0; ip<t_track.nplanes; ip++)cout<< t_track.GetSensor(ip)+1;
3434                      cout << " >> ADD ";
3435                  }
3436    
3437                  _caloTj->DoTrack(t_track.al,t_track.zini);
3438                  for(int ip=0; ip<_caloTj->npoint; ip++){
3439                      t_track.xGF[ip] = _caloTj->x[ip];
3440                      t_track.yGF[ip] = _caloTj->y[ip];
3441                  }
3442    
3443                  //              cout << "ADD "<<endl;
3444                  trackCandidates[mapIndex].insert(pair<int,ExtTrack>(t_track.GetNX()+t_track.GetNY(), ExtTrack(t_track)));
3445                  if(t_track.GetNX()+t_track.GetNY()>5 && t_track.chi2<chi2Min) chi2Min=t_track.chi2;
3446                  iCand++;
3447              }//end add-candidate condition
3448                
3449          }//end first loop on candidates
3450          if(_debug)cout <<endl<< "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
3451    
3452    
3453    
3454          //=================================================================
3455          // ----------------------------------------------------------------
3456          // chi2 selection
3457          // ----------------------------------------------------------------
3458          //=================================================================
3459          pair <multimap<int,ExtTrack>::iterator, multimap<int,ExtTrack>::iterator> nxny[NEXTVIEWS];
3460    //  for(int inxny=0; inxny<NEXTVIEWS; inxny++){
3461          bool NMTI =false;//not-minimal track inserted
3462          for(int inxny=NEXTVIEWS-1; inxny>=nClFixX+nClFixY; inxny--){
3463              nxny[inxny] = trackCandidates[mapIndex].equal_range(inxny);
3464    //    cout << " --- n "<<inxny<<endl;
3465              for ( multimap<int,ExtTrack>::iterator it = nxny[inxny].first; it!=nxny[inxny].second; ++it){
3466                  ExtTrack &trackC = (*it).second;
3467    //      cout << " NX+NY "<< trackC.GetNX()+trackC.GetNY() <<" chi2 "<<trackC.chi2<< " R(GV) "<<trackC.GetRigidity()<<endl;
3468                  //            if( trackC.chi2==0 )continue;
3469                  if( (trackC.GetNX()+trackC.GetNY() > 5 && trackC.chi2>chi2Min  ) ||
3470                      (trackC.GetNX()+trackC.GetNY() == 5 && NMTI ) ||
3471                      (trackC.GetNX()+trackC.GetNY() == 5 && trackC.chi2>1. ) ||
3472                      false)continue;      
3473                  
3474                  trackCandidates[mapIndex+1].insert(pair<int,ExtTrack>(inxny, ExtTrack(trackC)));
3475    //      cout << " insert "<<endl;
3476                  
3477                  if( trackC.GetNX()+trackC.GetNY() > 5 )NMTI=true;
3478              }
3479          }
3480          mapIndex++;
3481          
3482          //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
3483          if(
3484              trackCandidates[mapIndex].size()==0 ||
3485              (trackCandidates[mapIndex].size()>1 && nViewCal < _alg_nViewCal ) ||
3486              false){
3487              mapIndex++;
3488              continue;//increment one calo plane
3489          }
3490          //-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
3491    
3492          if(_debug){    
3493        
3494              cout<< " Selected-track SUMMARY:"<<endl;
3495    
3496    //        if(mapIndex>0)trackCandidates[mapIndex-1].clear();
3497              trkBegin = trackCandidates[mapIndex].begin();
3498              trkEnd   = trackCandidates[mapIndex].end();
3499              for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
3500                  ExtTrack &trackCand = (*trkIter).second;
3501                  cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
3502                  cout <<" | al: ";
3503                  for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
3504                  cout <<" | ";
3505                  cout << " X: ";
3506                  //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.XGood(ip) ? "o" : "-");
3507                  for(int ip=0; ip<trackCand.nplanes; ip++){
3508                      if(trackCand.GetClusterX_ID(ip)>=0)
3509                          cout << setw(2)<<trackCand.GetClusterX_ID(ip);
3510                      else
3511                          cout << "__";
3512                  }
3513                  cout << " | ";
3514                  cout << " Y: ";
3515                  //      for(int ip=0; ip<NEXTPLANES; ip++)cout << (trackCand.YGood(ip) ? "o" : "-");
3516                  for(int ip=0; ip<trackCand.nplanes; ip++){
3517                      if(trackCand.GetClusterY_ID(ip)>=0)
3518                          cout << setw(2)<<trackCand.GetClusterY_ID(ip);
3519                      else
3520                          cout << "__";
3521                  }
3522                  cout << " S:";for(int ip=0; ip<trackCand.nplanes; ip++)cout<< trackCand.GetSensor(ip)+1;
3523                  cout <<"  MDR="<<1./sqrt(trackCand.coval[4][4])<<endl;
3524                  cout << endl;
3525                  // for(int ip=0; ip<trackCand.nplanes; ip++)
3526                  //        if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";       cout << endl;
3527    
3528              }
3529          }
3530    
3531    
3532          //=================================================================
3533          // ----------------------------------------------------------------
3534          // save tracks
3535          // ----------------------------------------------------------------
3536          //=================================================================
3537    //      cout << mapIndex << endl;
3538          if(mapIndex>=maxSize)cout <<" ORRORE!!!! mapIndex"<<mapIndex<<" >= "<<maxSize<<endl;
3539    //      if(mapIndex>0)trackCandidates[mapIndex-1].clear();
3540          trkBegin = trackCandidates[mapIndex].begin();
3541          trkEnd   = trackCandidates[mapIndex].end();
3542          iCand = 0;
3543          // TClonesArray &tracks = *_trkArray;
3544          // tracks.Clear();
3545          for( multimap<int,ExtTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){        
3546              ExtTrack trackCand = ((*trkIter).second);
3547    //        cout << iCand << " "<< trackCand.al[4] <<endl;
3548              //    trackCand.seqno = iCand;
3549    //        trackCand.Dump();
3550              new(tracks[iCand]) ExtTrack(trackCand);
3551              iCand++;
3552          };
3553          break;
3554    
3555    
3556    
3557      }// end loop over n.calo planes to be included
3558    
3559    
3560    
3561    
3562    
3563    }
3564    
3565    
3566    /**
3567     *
3568     */
3569    void ExtTrkingAlg::ProcessEvent(Bool_t force){
3570        
3571        if(_whichAlg == 0)                           ProcessEvent0(force);
3572        else if(_whichAlg >= 100 && _whichAlg < 144) ProcessEvent1(force);
3573        else if(_whichAlg >= 200 && _whichAlg < 244) ProcessEvent2(force);
3574        else{
3575            cout << " Algorythm id "<< _whichAlg <<" not implemented "<<endl;
3576        }
3577    
3578    };
3579    /**
3580     *
3581     */
3582    void ExtTrkingAlg::Dump(){
3583    
3584      cout << "|--|--|--|--|--|--|--|--|--|--|--|"<<endl;
3585      cout << "ExtTrkingAlg::Dump() "<<endl<<endl;
3586      cout << "  algorythm "<<setw(10)<<_whichAlg <<endl;
3587      cout << "  debug     "<<setw(10)<< _debug <<endl;
3588      cout << "  TrkLevel1 "<<setw(10)<<_trk_l1 <<endl;
3589      cout << "  TrkLevel2 "<<setw(10)<<_trk_l2 <<endl;
3590      cout << "  CaloLevel1"<<setw(10)<<_cal_l1 <<endl;
3591      cout << "  CaloLevel2"<<setw(10)<<_cal_l2 <<endl;
3592      cout << "  ToFLevel2 "<<setw(10)<<_tof_l2 <<endl;
3593      cout << "Pre-selection parameters "<<endl;
3594      cout << "  nClstrMAX  "<<setw(10)<<_sel_nClstrMAX <<endl;
3595      cout << "  nPlaneXMIN "<<setw(10)<<_sel_nPlaneXMIN <<endl;
3596      cout << "  nPlaneYMIN "<<setw(10)<<_sel_nPlaneYMIN <<endl;
3597      cout << "Algorythm parameters "<<endl;
3598      cout << "  nClFixX    "<<setw(10)<<_alg_nClFixX <<endl;
3599      cout << "  nClFixY    "<<setw(10)<<_alg_nClFixY <<endl;
3600      cout << "  nTrackMAX  "<<setw(10)<<_alg_nTrackMAX <<endl;
3601      cout << "  nViewCal   "<<setw(10)<<_alg_nViewCal <<endl;
3602      cout << "|--|--|--|--|--|--|--|--|--|--|--|"<<endl;
3603      
3604    }
3605    
3606    
3607    
3608    ////////////////////////////////////////////////////////////////////////
3609    //
3610    //
3611    //
3612    //
3613    ////////////////////////////////////////////////////////////////////////
3614    //  int GetSensor(int view, int sis){
3615    //    if(sis<0||sis>8)return -1;
3616    //    if(view<0||view>1)return -1;
3617    //    if(!view)return (int)(sis/3);
3618    //    else return sis%3;
3619    //  };
3620    //  int GetLadder(int view, int sis){
3621    //    if(sis<0||sis>8)return -1;
3622    //    if(view<0||view>1)return -1;
3623    //    if(!view)return sis%3;  
3624    //    else return (int)(sis/3);
3625    //  };
3626    //  int GetSiSensor(int view, int l, int s){
3627    //    if(view<0||view>1)return -1;
3628    //    if(s<0||s>2)return -1;
3629    //    if(!view) return 3*s+l;
3630    //    else      return 3*l+s;
3631    //   };
3632    
3633      int CaloStripRoto::GetSensor(){
3634          return GetCaloSensor(GetView(),GetSiSensor());
3635      };
3636      int CaloStripRoto::GetLadder(){
3637          return GetCaloLadder(GetView(),GetSiSensor());
3638      };
3639    
3640     void CaloStripRoto::ResetAligParams(){
3641       fPitch = 1.;
3642       shift[0]=0.; //L
3643       shift[1]=0.; //S  
3644       shift[2]=0.; //Z
3645       alpha = 0.;
3646      
3647    //   siSensor = -1;;
3648    
3649     };
3650    /*1
3651     *
3652    
3653     *
3654     */
3655    bool CaloStripRoto::SensorContains(float x, float y){
3656    
3657        float toll = 0.2;//cm
3658    
3659    //    cout << "S -- SetStrip("<<GetLadder()*32+0<<")"<<endl;
3660        
3661    //    cout << " extst.SetStrip("<<GetLadder()*32+0<<") "<<endl;
3662        SetStrip(GetLadder()*32+0);
3663        float xM0,yM0;
3664        float d0 = GetDistanceTo(x,y,xM0,yM0);
3665    
3666        if( yM0-toll > fYB && yM0 > fYA )return false;
3667        if( yM0+toll < fYB && yM0 < fYA )return false;
3668    
3669        if(  xM0-toll > fXB && xM0 > fXA )return false;
3670        if(  xM0+toll < fXB && xM0 < fXA )return false;
3671    
3672    //    cout << "D -- SetStrip("<<GetLadder()*32+31<<")"<<endl;
3673    //    cout << " extst.SetStrip("<<GetLadder()*32+31<<") "<<endl;
3674        SetStrip(GetLadder()*32+31);
3675        float xM31,yM31;
3676        float d31 = GetDistanceTo(x,y,xM31,yM31);
3677    
3678    
3679        if( d0+d31 > 31*fPitch*0.244+2.*toll)return false;
3680    
3681        return true;
3682    
3683    };
3684    
3685     CaloStripRoto::CaloStripRoto(int view, int plane, int sisensor, Bool_t usemechanicalalignement){
3686        
3687         SetView(view);
3688         SetPlane(plane);
3689         SetSiSensor(sisensor);
3690    
3691         st = CaloStrip(usemechanicalalignement);
3692         ResetAligParams();
3693         if(!usemechanicalalignement){
3694    //       cout <<" SETTA I PARAMETRI DI ALLINEAMNETO  -- da implementare"<<endl;
3695         }
3696        
3697    //     Set(view,plane,0,GetSiSensor());
3698     } ///< Default Constructor.
3699    
3700    
3701    
3702     //         relative to pamela reference system (convenzione emiliano):
3703     //
3704     //        y ^
3705     //          |
3706     //          | 6 7 8 sisensor
3707     //          | 3 4 5  
3708     //          | 0 1 2
3709     //          ------------> x
3710     //
3711     //        relative to silicon sensor          
3712     //
3713     //        S ^
3714     //          |   |
3715     //        2 |   |  <- strip direction
3716     // sensor 1 |   |
3717     //        0 |...|...
3718     //          ------------> L
3719     //            0 1 2
3720     //            ladder
3721     //
3722     //
3723    void CaloStripRoto::Set(Int_t view, Int_t plane, float strip, int sisensor){
3724    
3725    
3726        SetView(view);
3727        SetPlane(plane);
3728    
3729    
3730        float size = 8.05; //(cm) larghezza del sensore + zone morte
3731      
3732      
3733       //----------------------------
3734       // Z-coordinate
3735       //----------------------------
3736    
3737        float coordCM_strip;
3738        if( strip >0. && strip < 95.){
3739            st.Set( view,plane,(int)strip);// istrip
3740            coordCM_strip = ( view ? st.GetY() : st.GetX() );
3741            coordCM_strip *= (1-strip+(int)strip);
3742            st.Set( view,plane,(int)strip+1);// istrip+1
3743            coordCM_strip += (strip-(int)strip)*( view ? st.GetY() : st.GetX() );
3744        }else{
3745    //      cout << "CaloStripRoto::Set( "<<view<<","<<plane<<","<<strip<<","<<sisensor<<")"<<endl;
3746    //      cout << " ---> ERROR: 0 =< strip =< 95 " << endl;
3747            if(strip>=95.)strip=95.;
3748            if(strip<=0.)strip=0.;
3749            st.Set( view,plane,(int)strip);// istrip
3750            coordCM_strip = ( view ? st.GetY() : st.GetX() );      
3751        }
3752        float coordCm_Z0 = st.GetZ();
3753        //----------------------------
3754        // strip-segment coordinates
3755        //----------------------------
3756        float coordCm_LA = coordCM_strip;
3757        float coordCm_LB = coordCM_strip;
3758        float coordCm_SA = (-1)*size*1.5; //full calorimeter size
3759        float coordCm_SB = (+1)*size*1.5; //full calorimeter size
3760        float coordCm_ZA = coordCm_Z0;
3761        float coordCm_ZB = coordCm_Z0;
3762    
3763      //----------------------------
3764      // rototraslation
3765      //----------------------------
3766       //if sisensor is assigned (0,...8) the strip can be rototraslated
3767       bool ROTO = (sisensor>=0&&sisensor<9);
3768    
3769       if(ROTO){
3770    
3771    
3772           if(        
3773               GetCaloSiSensor(view,(int)(strip/32),0) != sisensor &&
3774               GetCaloSiSensor(view,(int)(strip/32),1) != sisensor &&
3775               GetCaloSiSensor(view,(int)(strip/32),2) != sisensor &&
3776               true){
3777               cout << " Strip "<<strip<<" does not belong to sensor "<<sisensor<<endl;
3778               cout << " input data are inconsistent --- output might be wrong "<<endl;    
3779           } ;
3780    
3781    
3782         int sensor = GetCaloSensor(view, sisensor);//0-1-2
3783         int ladder = GetCaloLadder(view, sisensor);//0-1-2
3784         if( ladder != (int)(strip / 32) ){
3785           cout << " void CaloStripRoto::Set(..) --- ERROR --- strip "<<strip<<" does not match Si-sensor "<<sisensor<<endl;
3786           return;
3787         }
3788    
3789         //-----------------------------
3790         // sensor origin coordinates
3791         //-----------------------------
3792        
3793         st.Set(         view, plane, ladder*32); //strip = 0-32-64
3794         float coordCm_L0 = ( view ? st.GetY() : st.GetX());  
3795         st.Set( (int)(!view), plane, sensor*32 );//strip = 0-32-64
3796         float coordCm_S0 = ( view ? st.GetX() : st.GetY());
3797    
3798         //set sensor size
3799         coordCm_SA = (-1)*size*0.5 + (sensor-1)*size;
3800         coordCm_SB = (+1)*size*0.5 + (sensor-1)*size;
3801    
3802         //coordinate relative al sensore
3803         coordCm_LA = coordCm_LA - coordCm_L0;
3804         coordCm_SA = coordCm_SA - coordCm_S0;
3805         coordCm_ZA = coordCm_ZA - coordCm_Z0;
3806         coordCm_LB = coordCm_LB - coordCm_L0;
3807         coordCm_SB = coordCm_SB - coordCm_S0;
3808         coordCm_ZB = coordCm_ZB - coordCm_Z0;
3809    
3810         //variazione di pitch
3811         coordCm_LA = coordCm_LA * fPitch;
3812         coordCm_LB = coordCm_LB * fPitch;
3813    
3814         //rotazione
3815         float LA = coordCm_LA - alpha * coordCm_SA;
3816         float SA = coordCm_SA + alpha * coordCm_LA;
3817         float LB = coordCm_LB - alpha * coordCm_SB;
3818         float SB = coordCm_SB + alpha * coordCm_LB;
3819    
3820         //traslazione
3821         coordCm_LA = LA         + (!view ? shift[0] : shift[1]);
3822         coordCm_SA = SA         + ( view ? shift[1] : shift[0]);
3823         coordCm_ZA = coordCm_ZA + shift[2];
3824         coordCm_LB = LB         + shift[0];
3825         coordCm_SB = SB         + shift[1];
3826         coordCm_ZB = coordCm_ZB + shift[2];
3827    
3828         //back to sistema di riferimento generale
3829        
3830         coordCm_LA = coordCm_LA + coordCm_L0;
3831         coordCm_SA = coordCm_SA + coordCm_S0;
3832         coordCm_ZA = coordCm_ZA + coordCm_Z0;
3833         coordCm_LB = coordCm_LB + coordCm_L0;
3834         coordCm_SB = coordCm_SB + coordCm_S0;
3835         coordCm_ZB = coordCm_ZB + coordCm_Z0;
3836    
3837       }
3838      
3839       fXA = (!view ? coordCm_LA : coordCm_SA);
3840       fYA = ( view ? coordCm_LA : coordCm_SA);
3841       fZA = coordCm_ZA;
3842    
3843       fXB = (!view ? coordCm_LB : coordCm_SB);
3844       fYB = ( view ? coordCm_LB : coordCm_SB);
3845       fZB = coordCm_ZB;
3846    
3847    
3848       st.Set(view,plane,(int)strip);// istrip
3849    
3850     };
3851     //
3852     //
3853     //
3854     //
3855     //
3856     //
3857     float CaloStripRoto::GetSpatialResolution(float def, float degx, float degy,  float beta){
3858      
3859       float res_dig = 0.24/sqrt(12.);
3860    
3861       if(def==0.)return res_dig;
3862      
3863       float res_ms = 0.;
3864      
3865       // Int_t GetView(){return (fView-1);} ///< Get strip view [0-1]
3866       // Int_t GetPlane(){return (fPlane-1);} ///< Get strip plane [0-21]
3867       // Int_t GetStrip(){return (fStrip-1);} ///< Get strip number [0-95]
3868      
3869      
3870       float rig = 1/def; if(rig<0)rig = -rig;
3871       float factor = 0.;
3872       factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*degx/180.),2);
3873       factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*degy/180.),2);
3874       factor += 1.;
3875       factor = TMath::Sqrt(factor);
3876      
3877       int view = st.GetPlane()*2 + 1 - st.GetView()  ;//cluster.view; //0... 43
3878      
3879       //  int nW = (int)((view + 1)/2);
3880       float dW = 0.74*factor;//X0
3881       float dW_cm = 0.26*factor;//cm
3882       float dSi_cm = 0.; //cm
3883       if(view && st.GetPlane()>0) dSi_cm = ( st.GetPlane()%2 ? 0.228 : 0.428);
3884       dSi_cm *= factor;
3885       // multiple scattering angle
3886       float theta = 0.;//ms
3887       theta = 0.0136/rig;
3888       theta = theta * sqrt(dW) * (1.+0.038*log(dW));
3889       float Vt = theta*theta;            //cov(theta,theta)
3890       float Vy = Vt * dW_cm * dW_cm /3.; //cov(y,y)
3891       float Vty = 0.87 * Vt * Vy;        //cov(theta,y)
3892       float Vc = Vy + dSi_cm * dSi_cm * Vt + 2. * dSi_cm * Vty;
3893       float rms = sqrt( Vc );        
3894       // if(isY)t_track.resy[ip] = sqrt(rms*rms + t_track.resy[ip]*t_track.resy[ip]);
3895       // else   t_track.resx[ip] = sqrt(rms*rms + t_track.resx[ip]*t_track.resx[ip]);
3896       // if(isY)          VMS[nW] = rms*rms;
3897       // if(_debug)cout <<endl<<view<<" nW "<<nW<<" res(MS) "<<rms<<" factor "<<factor;
3898      
3899       if(view>=3){
3900         cout << " --> NB resolution not accurate for layer "<<view<<endl;;
3901         cout <<"st.GetPlane() "<<st.GetPlane()<<endl;
3902         cout <<"st.GetView() "<<st.GetView()<<endl;
3903       }
3904       res_ms = rms;
3905      
3906       return sqrt(res_ms*res_ms + res_dig*res_dig);
3907        
3908       };
3909    /**
3910     *
3911     *
3912     *
3913     *
3914     */
3915    float CaloStripRoto::GetDistanceTo(float xP, float yP, float& xM, float &yM ){
3916    
3917    
3918        bool isX = fabs(fXA-fXB) < fabs(fYA-fYB);
3919    
3920    
3921        if(isX){
3922            float beta = (fXB-fXA)/(fYB-fYA);
3923            float alpha = fXA - beta * fYA;
3924            yM = ( yP + beta*xP - beta*alpha )/(1+beta*beta);
3925            xM = alpha + beta * yM;
3926        }else{
3927            float beta = (fYB-fYA)/(fXB-fXA);
3928            float alpha = fYA - beta * fXA;
3929            xM = ( xP + beta*yP - beta*alpha )/(1+beta*beta);
3930            yM = alpha + beta * xM;
3931        }
3932        
3933        return sqrt((xM-xP)*(xM-xP) + (yM-yP)*(yM-yP));
3934    
3935    };

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4

  ViewVC Help
Powered by ViewVC 1.1.23