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

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

  ViewVC Help
Powered by ViewVC 1.1.23