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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (hide annotations) (download)
Thu Jul 24 12:39:50 2014 UTC (10 years, 4 months ago) by pam-ts
Branch: MAIN
Changes since 1.4: +83 -26 lines
TrkCore bug fixed

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

  ViewVC Help
Powered by ViewVC 1.1.23