/[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.1 - (hide annotations) (download)
Thu Feb 27 11:24:43 2014 UTC (10 years, 9 months ago) by pam-fi
Branch: MAIN
Added new tracking algorythm

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     */
13     ExtTrkingAlg::ExtTrkingAlg(Int_t id){
14    
15     _whichAlg = id;
16    
17     SetDebug(false);
18    
19     double pars[] = {20.,3.,2.};
20     SetSelectionParams(pars);
21    
22    
23     double para[] = {3.,2.,400.};
24     SetAlgorythmParams(para);
25    
26     if(id == 0){
27    
28     cout << "ExtTrkingAlg::ExtTrkingAlg("<<id<<")"<<endl;
29     cout << "Creating array of TrkTrack objects "<<endl;
30     _trkArray = new TClonesArray("TrkTrack");
31    
32     }else{
33    
34     cout << "ExtTrkingAlg(Int_t id) -- algorythm id="<<id<<" not valid "<<endl;
35     cout << "--> Track array not created "<<endl;
36    
37     }
38    
39     Clear();
40    
41     };
42    
43     /**
44     * Reset the algorythm flags and output
45     */
46     void ExtTrkingAlg::Reset(){
47    
48     if(_trkArray)_trkArray->Clear("C");
49    
50    
51    
52     };
53    
54     /**
55     * Clear the event and reset the algorythm
56     */
57     void ExtTrkingAlg::Clear(Option_t* option){
58    
59     SetTrkLevel1();
60     SetTrkLevel2();
61    
62     SetToFLevel2();
63    
64     SetCaloLevel1();
65     SetCaloLevel2();
66    
67     Reset();
68    
69     };
70     /**
71     * Set selection parameters
72     */
73     void ExtTrkingAlg::SetSelectionParams(double* par){
74    
75     _sel_nClstrMAX = (Int_t)par[0]; // selection parameter: maximum number of cluster
76     _sel_nPlaneXMIN = (Int_t)par[1]; // selection parameter: minimum number of hit x-views
77     _sel_nPlaneYMIN = (Int_t)par[2]; // selection parameter: minimum number of hit y-views
78    
79     };
80     /**
81     * Set algorythm parameters
82     */
83     void ExtTrkingAlg::SetAlgorythmParams(double* par){
84    
85     _alg_nClFixX = (Int_t)par[0]; // algorythm parameter: n.hits required on X view
86     _alg_nClFixY = (Int_t)par[1]; // algorythm parameter:n.hits required on X view
87     _alg_nTrackMAX = (Int_t)par[2]; // algorythm parameter: maximum num. of track candidates
88    
89     };
90    
91     /**
92     * Get the track array
93     */
94     TClonesArray* ExtTrkingAlg::GetTrackArray(Bool_t reset){
95    
96     if( reset ) Reset();
97     // if( !_procDone )ProcessEvent();
98    
99     return _trkArray;
100    
101     };
102     /**
103     * Pre-select the event
104     */
105     Bool_t ExtTrkingAlg::CheckEvent(){
106    
107     if(!_trk_l2)return true;
108    
109     TrkLevel2 *trkl2 = _trk_l2;//event->GetTrkLevel2();
110     int nClstrMAX = _sel_nClstrMAX; //maximum number of cluster
111     int nPlaneXMIN = _sel_nPlaneXMIN;
112     int nPlaneYMIN = _sel_nPlaneYMIN;
113     /////////////////////////////////////////////////////////////////////
114     /// count number of hit planes
115     /////////////////////////////////////////////////////////////////////
116     int nHitX[] = {0,0,0,0,0,0};
117     int nHitY[] = {0,0,0,0,0,0};
118     for(int ix=0; ix<trkl2->nclsx(); ix++)nHitX[trkl2->GetSingletX(ix)->plane-1]++;
119     for(int iy=0; iy<trkl2->nclsy(); iy++)nHitY[trkl2->GetSingletY(iy)->plane-1]++;
120     /////////////////////////////////////////////////////////////////////
121     /// set minimum condition 3x+2y
122     /////////////////////////////////////////////////////////////////////
123     int nPlaneX=0;
124     for(int ip=0; ip<6; ip++)if(nHitX[ip])nPlaneX++;
125     int nPlaneY=0;
126     for(int ip=0; ip<6; ip++)if(nHitY[ip])nPlaneY++;
127     /////////////////////////////////////////////////////////////////////
128     /// dump selection
129     /////////////////////////////////////////////////////////////////////
130     if(_debug){
131     cout << " n.tracks "<<trkl2->GetNTracks()<<endl;
132     cout << " n.x-singles "<<trkl2->nclsx()<<endl;
133     cout << " n.y-singles "<<trkl2->nclsy()<<endl;
134     cout << " n.x-planes with singles "<<nPlaneX<<endl;
135     cout << " n.y-planes with singles "<<nPlaneY<<endl;
136     }
137     ////////////////////////////////////////////////////////////////////////
138     //
139     // check condition for retracking
140     //
141     ////////////////////////////////////////////////////////////////////////
142     if(trkl2->GetNTracks()==0 && nPlaneX<nPlaneXMIN) return false;
143     if(trkl2->GetNTracks()==0 && nPlaneY<nPlaneYMIN) return false;
144     if(trkl2->nclsy()+trkl2->nclsx() > nClstrMAX) return false;
145    
146     return true;
147    
148     };
149    
150    
151     /**
152     * Process the event
153     */
154     void ExtTrkingAlg::ProcessEvent(Bool_t force){
155    
156     if(_debug)cout << " |---------------------------------------------------| "<<endl;
157     if(_debug)cout << " void ExtTrkingAlg::ProcessEvent("<<force<<") "<<endl;
158    
159    
160     if(_debug && !_trk_l1)cout << " Missing TrkLevel1 object --- ciao! "<<endl;
161     if(!_trk_l1)return;
162    
163    
164     //===========================
165     // pre-selection
166     //===========================
167    
168     Bool_t RETRACK = CheckEvent();
169    
170     if(_debug)cout <<" Pre-selection result: "<<RETRACK<<endl;
171     if( !force && !RETRACK )return;
172    
173    
174     //===========================
175     // tracking
176     //===========================
177    
178     TClonesArray &tracks = *_trkArray;
179     tracks.Clear();
180     multimap<int,TrkTrack> trackCandidates[12]; // map: last filled view vs TrkTrack
181     int mapIndex = 0;
182     multimap<int,TrkTrack>::iterator trkBegin;
183     multimap<int,TrkTrack>::iterator trkEnd;
184     int iCand;
185    
186     // ------------------------------------------------
187     // fill a map VIEW(0-11)-CLUSTER
188     // ------------------------------------------------
189     multimap<int,int> clusterMap;
190     if(_debug)cout<<"-------- CLUSTERs --------"<<endl;
191     for(int icl=0; icl<_trk_l1->nclstr(); icl++){
192     clusterMap.insert( pair<int,int>( _trk_l1->GetCluster(icl)->view-1, icl ));
193     if(_debug)cout << icl << " - view "<<_trk_l1->GetCluster(icl)->view<<" ladder "<<_trk_l1->GetCluster(icl)->GetLadder()<<" strip "<<_trk_l1->GetCluster(icl)->maxs<<endl;
194     }
195     // ------------------------------------------------
196     // define iterators over clusters on the same view
197     // ------------------------------------------------
198     pair <multimap<int,int>::iterator, multimap<int,int>::iterator> ret[12];
199     multimap<int,int>::iterator it[12];
200     for(int iv=0; iv<12; iv++)ret[iv] = clusterMap.equal_range(iv);
201    
202     //====================================================================
203     // PARAMETERS
204     //====================================================================
205     // ------------------------------------------------
206     // number of points per candidate
207     // ------------------------------------------------
208     int nClFixX = _alg_nClFixX; // n.hits required on X view
209     int nClFixY = _alg_nClFixY; // n.hits required on X view
210     int nTrackMAX = _alg_nTrackMAX; //
211     //==================================================
212     // ------------------------------------------------
213     // start: one candidate for each cluster
214     // ------------------------------------------------
215     //==================================================
216     if(_debug)cout<<"-------- MINIMAL --------"<<endl;
217     for(int iv=0; iv<12; iv++){ //loop over all views
218    
219     if( !(iv%2) && (int)(0.5*iv) > 6-nClFixX)continue;
220     if( (iv%2) && (int)(0.5*iv) > 6-nClFixY)continue;
221    
222     for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view
223     for (int is=0; is<2; is++){ //loop over sensors
224     int icl = it[iv]->second;
225     TrkCluster *cluster = _trk_l1->GetCluster(icl);
226     int ladder = cluster->GetLadder()-1;
227     int plane = cluster->GetPlane()-1;
228     bool isY = (cluster->view)%2;
229     bool isX = !isY;
230     //
231     // if(_debug && isX )cout << " ^^^^^^^ X plane "<<plane<<" ladder "<<ladder<<" strip "<<cluster->maxs <<endl;
232     // if(_debug && isY )cout << " ^^^^^^^ Y plane "<<plane<<" ladder "<<ladder<<" strip "<<cluster->maxs <<endl;
233    
234     TrkTrack t_track = TrkTrack();
235    
236     if( isX )t_track.SetXGood(plane,icl+1,ladder,is);
237     if( isY )t_track.SetYGood(plane,icl+1,ladder,is);
238    
239     trackCandidates[mapIndex].insert(pair<int,TrkTrack>(iv, TrkTrack(t_track)));
240    
241     }
242     }
243     };
244    
245     //==============================================================
246     // -------------------------------------------------------------
247     // next: add views to each candidate, up to nClFixX+nClFixY
248     // -------------------------------------------------------------
249     //==============================================================
250     for( ;mapIndex < nClFixX+nClFixY-1 ;++mapIndex){
251    
252     trkBegin = trackCandidates[mapIndex].begin();
253     trkEnd = trackCandidates[mapIndex].end();
254     if(_debug)cout<<"MAP "<<mapIndex<<" -----------"<<endl;
255    
256     for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
257    
258     TrkTrack trackCand = (*trkIter).second;
259     int lastView = (*trkIter).first; //last inserted view
260     int lastPlane = (int)(0.5*lastView);
261     int lastLadder = trackCand.GetLadder(lastPlane);
262     int lastSensor = trackCand.GetSensor(lastPlane);
263    
264     // if(_debug)cout<<"MAP "<<mapIndex<<" last v "<<lastView<<" NX "<<trackCand.GetNX()<<" NY "<<trackCand.GetNY()<<endl;
265    
266     for(int iv=lastView+1; iv<12; iv++){ //loop over next views
267    
268     if( (iv%2) && (int)(0.5*iv) > 6-nClFixX+trackCand.GetNX())continue; //not enough x-view left
269     if( !(iv%2) && (int)(0.5*iv) > 6-nClFixY+trackCand.GetNY())continue; //not enough y-view left
270    
271     if( (iv%2) && trackCand.GetNX()==nClFixX )continue;//ok?
272     if( !(iv%2) && trackCand.GetNY()==nClFixY )continue;//ok?
273    
274     for ( it[iv]=ret[iv].first; it[iv]!=ret[iv].second; ++it[iv]){ //loop over clusters on the view
275    
276     int icl = it[iv]->second;
277     TrkCluster *cluster = _trk_l1->GetCluster(icl);
278     int ladder = cluster->GetLadder()-1;//(_trk_l1->GetCluster(icl)->maxs-1)%1024 + 1;
279     int plane = cluster->GetPlane()-1;
280     bool isY = (cluster->view)%2;
281     bool isX = !isY;
282    
283     if( plane==lastPlane && ladder!=lastLadder)continue;
284    
285    
286     for (int is=0; is<2; is++){ //loop over sensors
287    
288     if( plane==lastPlane && is!=lastSensor)continue;
289    
290     // if(_debug)cout<<"MAP "<<mapIndex<<" last v "<<lastView;
291     // if(_debug)cout<<" +cl "<<icl<<" s"<<is<<endl;
292    
293     TrkTrack t_track = TrkTrack();
294     trackCand.Copy(t_track); //copy the candidate parameters...
295    
296     if( isX )t_track.SetXGood(plane,icl+1,ladder,is); //...add a point...
297     if( isY )t_track.SetYGood(plane,icl+1,ladder,is);
298    
299     trackCandidates[mapIndex+1].insert(pair<int,TrkTrack>(iv, TrkTrack(t_track))); //...and store a new candidate
300     };//end loop over sensors
301     };//end loop over clusters on the view iv
302     };//endl loop over views
303     }//end loop over candidates
304     }//end iterations
305    
306    
307     if( trackCandidates[mapIndex].size() > nTrackMAX ){
308     if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<" > "<<nTrackMAX<<endl;
309     return;//to many candidates
310     }
311     if(_debug){
312     trkBegin = trackCandidates[mapIndex].begin();
313     trkEnd = trackCandidates[mapIndex].end();
314     for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
315     TrkTrack &trackCand = (*trkIter).second;
316     cout << " >> ";
317     cout << " X: ";
318     for(int ip=0; ip<6; ip++)
319     if(trackCand.GetClusterX_ID(ip)>=0)cout << trackCand.GetClusterX_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";
320     cout << " Y: ";
321     for(int ip=0; ip<6; ip++)
322     if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";
323     cout << endl;
324     }
325     }
326     if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
327    
328     //==============================================================
329     // -------------------------------------------------------------
330     // first selection of track candidates, based on chi2
331     // -------------------------------------------------------------
332     //==============================================================
333     if(_debug)cout<<"-------- CHI2 --------"<<endl;
334     trkBegin = trackCandidates[mapIndex].begin();
335     trkEnd = trackCandidates[mapIndex].end();
336     mapIndex++;
337     iCand = 0;
338     for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
339     int lastView = (*trkIter).first;
340     TrkTrack &trackCand = (*trkIter).second;
341     int ifail=0;
342     trackCand.FitReset();
343     trackCand.Fit(0.,ifail,0,1);//re-evaluate cluster position from level1
344     if(ifail!=0)trackCand.FitReset();
345    
346     // -----------------------------------
347     // first selection of track candidates
348     // -----------------------------------
349     if(TMath::IsNaN(trackCand.chi2))continue;
350     if(trackCand.chi2 <= 0. || trackCand.chi2 > 1.e7)continue;
351     // if(trackCand.chi2 < 0. || trackCand.chi2 > 10. )continue;
352     iCand++;
353    
354     if(_debug){
355     cout<< " TRK chi2 "<<setw(13)<<trackCand.chi2;
356     cout <<" | al: ";
357     for(int i=0; i<5; i++)cout <<setw(10)<< trackCand.al[i];
358     cout <<" | ";
359     for(int ip=0; ip<6; ip++)
360     if(trackCand.GetClusterX_ID(ip)>=0)cout << trackCand.GetClusterX_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|";
361     cout << " Y: ";
362     for(int ip=0; ip<6; ip++)
363     if(trackCand.GetClusterY_ID(ip)>=0)cout << trackCand.GetClusterY_ID(ip)<<":"<<trackCand.GetSensor(ip)<<"|"; cout << endl;
364     }
365    
366     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
367     // evaluated coordinates (to define GF)
368     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
369     int ngf = TrkParams::nGF;
370     float *zgf = TrkParams::zGF;
371     Trajectory tgf = Trajectory(ngf,zgf);
372     tgf.DoTrack(trackCand.al);//<<<< integrate the trajectory
373     for(int ip=0; ip<ngf; ip++){
374     trackCand.xGF[ip] = tgf.x[ip];
375     trackCand.yGF[ip] = tgf.y[ip];
376     }
377    
378     trackCandidates[mapIndex].insert(pair<int,TrkTrack>(lastView, TrkTrack(trackCand)));
379    
380     };
381     if(trackCandidates[mapIndex].size()==0)return;//no good candidates
382    
383     if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
384     //=================================================================
385     // ----------------------------------------------------------------
386     // CALORIMETER
387     // ----------------------------------------------------------------
388     //=================================================================
389     if(_cal_l1){
390     if(_debug)cout<<"-------- CALORIMETER --------"<<endl;
391     // -----------------------------------
392     // evaluate calorimeter variables
393     // -----------------------------------
394     float caloCoord[2][22][96];
395     float caloZ[44];
396     for(int view=0; view<2; view++){
397     for(int plane=0; plane<22; plane++){
398     for(int strip=0; strip<96; strip++){
399     CaloStrip st = CaloStrip(false);
400     st.Set(view,plane,strip);
401     caloCoord[view][plane][strip]=(view ? st.GetY() : st.GetX());
402     caloZ[plane*2+(view?0:1)]=st.GetZ();
403     }
404     }
405     }
406     Trajectory caloTj = Trajectory(44,caloZ);
407     // for(int iz=0; iz<44;iz++)cout<<caloZ[iz]<<endl;
408     //test calo
409     // float caloEvent[2][22][96];
410     CaloLevel1* l1=_cal_l1;//event->GetCaloLevel1();
411     // if(!l1)return;
412     vector<float> caloChi2;
413     vector<int> caloHit; int caloHitMAX=0;
414     vector<float> caloMip;
415     trkBegin = trackCandidates[mapIndex].begin();
416     trkEnd = trackCandidates[mapIndex].end();
417     for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
418     TrkTrack trackCand = ((*trkIter).second);
419    
420     trackCand.DoTrack(&caloTj);//integrate the trajectory through the calorimeter
421    
422     float chi2 = 0.;
423     int nhit = 0;
424     float cmip = 0.;
425     for(Int_t ih=0; ih<l1->istrip; ih++){
426     Int_t view = -1;
427     Int_t plane = -1;
428     Int_t strip = -1;
429     Float_t mip = l1->DecodeEstrip(ih,view,plane,strip);
430     if(view<0)continue;
431     if(plane<0)continue;
432     if(strip<0)continue;
433     if(mip<=0)continue;
434     float dxy = caloCoord[view][plane][strip]-(view?caloTj.y[plane*2+(view?0:1)]:caloTj.x[plane*2+(view?0:1)]);
435     if(plane<11) chi2+= mip * TMath::Power(dxy,2.);
436     // caloZ[plane*2+(view?0:1)]=st.GetZ();
437     if( fabs(dxy) < 2. && plane<11)nhit++;
438     if( fabs(dxy) < 2. && plane<11)cmip+=mip;
439     };
440     if(l1->istrip>0)chi2 = chi2/l1->istrip;
441    
442     caloHit.push_back(nhit);
443     if(nhit>caloHitMAX)caloHitMAX=nhit;
444    
445     if(_debug){
446     cout << " CALO ";
447     // cout << " chi2 "<< chi2 <<" ";
448     cout << " nhit (< 11th plane) "<< nhit <<" ";
449     // cout << " mip "<< cmip <<" ";
450     cout <<endl;
451     };
452     }
453    
454     //=================================================================
455     // ----------------------------------------------------------------
456     // selection of track candidates, based on n.calorimeter hits
457     // ----------------------------------------------------------------
458     //=================================================================
459     trkBegin = trackCandidates[mapIndex].begin();
460     trkEnd = trackCandidates[mapIndex].end();
461     mapIndex++;
462     iCand = 0;
463     for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
464     int lastView = (*trkIter).first;
465     TrkTrack &trackCand = (*trkIter).second;
466     if(caloHit[iCand]==caloHitMAX)trackCandidates[mapIndex].insert(pair<int,TrkTrack>(lastView, TrkTrack(trackCand)));
467     iCand++;
468     }
469     if(trackCandidates[mapIndex].size()==0)return;//no good candidates
470     if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
471    
472     }
473    
474     //=================================================================
475     // ----------------------------------------------------------------
476     // TOF
477     // ----------------------------------------------------------------
478     //=================================================================
479     if(_tof_l2){
480     if(_debug)cout<<"-------- TOF --------"<<endl;
481     trkBegin = trackCandidates[mapIndex].begin();
482     trkEnd = trackCandidates[mapIndex].end();
483     ToFLevel2* tl2=_tof_l2;//event->GetToFLevel2();
484     // if(!tl2)return;
485     iCand =0;
486     vector<int> tofHit; int tofHitMAX=0;
487     for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
488     // int lastView = (*trkIter).first;
489     TrkTrack &trackCand = (*trkIter).second;
490    
491     int nhit =0;
492     for(int iToF=0; iToF<6; iToF++){
493    
494     int iGF = iToF;
495     if(iToF>3)iGF = iToF+8;
496     int iPaddle = tl2->GetPaddleIdOfTrack( trackCand.xGF[iGF], trackCand.yGF[iGF], iToF);
497     if( tl2->HitPaddle(iToF,iPaddle) )nhit++;
498    
499     }
500     iCand++;
501    
502     tofHit.push_back(nhit);
503     if(nhit>tofHitMAX)tofHitMAX=nhit;
504    
505     if(_debug){
506     cout << " TOF nhit "<< nhit <<" ";
507     cout <<endl;
508     }
509    
510     }
511     //=================================================================
512     // ----------------------------------------------------------------
513     // selection of track candidates, based on n.tof hits
514     // ----------------------------------------------------------------
515     //=================================================================
516     trkBegin = trackCandidates[mapIndex].begin();
517     trkEnd = trackCandidates[mapIndex].end();
518     mapIndex++;
519     iCand = 0;
520     for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
521     int lastView = (*trkIter).first;
522     TrkTrack &trackCand = (*trkIter).second;
523     if(tofHit[iCand]==tofHitMAX)trackCandidates[mapIndex].insert(pair<int,TrkTrack>(lastView, TrkTrack(trackCand)));
524     iCand++;
525     }
526     if(trackCandidates[mapIndex].size()==0)return;//no good candidates
527     if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
528    
529     }
530    
531    
532     //=================================================================
533     // ----------------------------------------------------------------
534     // TRACK MERGING
535     // ----------------------------------------------------------------
536     //=================================================================
537     if(_debug)cout<<"-------- MERGING --------"<<endl;
538     trkBegin = trackCandidates[mapIndex].begin();
539     trkEnd = trackCandidates[mapIndex].end();
540     mapIndex++;
541     float chi2Min=1e10;
542     iCand = 0;
543     for( multimap<int,TrkTrack>::iterator trkIter1 = trkBegin; trkIter1!=trkEnd; ++trkIter1){ //loop over tracks (1)
544    
545     TrkTrack &trackCand1 = (*trkIter1).second;// get candidate 1...
546    
547     TrkTrack t_track = TrkTrack();// ...create a new TrkTrack instance...
548     trackCand1.Copy(t_track); // ...and make a copy
549    
550     bool HBM = false; //has been merged
551    
552     for( multimap<int,TrkTrack>::iterator trkIter2 = trkBegin; trkIter2!=trkEnd; ++trkIter2){//loop over tracks (2)
553    
554     TrkTrack &trackCand2 = (*trkIter2).second; // get candidate 2...
555    
556     // if(_debug){
557     // cout <<" X: ";
558     // for(int ip=0; ip<6; ip++)
559     // if(t_track.GetClusterX_ID(ip)>=0)cout << t_track.GetClusterX_ID(ip)<<":"<<t_track.GetSensor(ip)<<"|";
560     // cout << " Y: ";
561     // for(int ip=0; ip<6; ip++)
562     // if(t_track.GetClusterY_ID(ip)>=0)cout << t_track.GetClusterY_ID(ip)<<":"<<t_track.GetSensor(ip)<<"|";
563     // cout << " + ";
564     // cout <<" X: ";
565     // for(int ip=0; ip<6; ip++)
566     // if(trackCand2.GetClusterX_ID(ip)>=0)cout << trackCand2.GetClusterX_ID(ip)<<":"<<trackCand2.GetSensor(ip)<<"|";
567     // cout << " Y: ";
568     // for(int ip=0; ip<6; ip++)
569     // if(trackCand2.GetClusterY_ID(ip)>=0)cout << trackCand2.GetClusterY_ID(ip)<<":"<<trackCand2.GetSensor(ip)<<"|";
570     // cout << " -----> ";
571     // }
572    
573     bool CBM = true;
574     for(int ip=0; ip<6; ip++){ //loop over planes
575     if( t_track.XGood(ip) &&
576     trackCand2.XGood(ip) &&
577     ( t_track.GetClusterX_ID(ip)!=trackCand2.GetClusterX_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
578     true )CBM=false; // if both have a x-hit and it is not the same ---> they cannot be merged
579     if(!CBM)break;
580     if( t_track.YGood(ip) &&
581     trackCand2.YGood(ip) &&
582     ( t_track.GetClusterY_ID(ip)!=trackCand2.GetClusterY_ID(ip) || t_track.GetSensor(ip)!=trackCand2.GetSensor(ip) ) &&
583     true )CBM=false; // if both have a y-hit and it is not the same ---> they cannot be merged
584     if(!CBM)break;
585     }
586     // if(_debug)cout << " can be merged ? "<<CBM<<endl;
587     if(!CBM)continue; //go to the next candidate
588     ///////////
589     //merging//
590     ///////////
591     for(int ip=0; ip<6; ip++){
592     if( !(t_track.XGood(ip)) && trackCand2.XGood(ip) )
593     t_track.SetXGood(ip,trackCand2.GetClusterX_ID(ip)+1,trackCand2.GetLadder(ip),trackCand2.GetSensor(ip));
594     if( !(t_track.YGood(ip)) && trackCand2.YGood(ip) )
595     t_track.SetYGood(ip,trackCand2.GetClusterY_ID(ip)+1,trackCand2.GetLadder(ip),trackCand2.GetSensor(ip));
596     };
597     HBM = true;
598     }//end loop over tracks (2)
599    
600    
601     // cout << "<< ";
602     // for(int ip=0; ip<6; ip++) cout <<setw(4)<< t_track.GetClusterX_ID(ip);
603     // for(int ip=0; ip<6; ip++) cout <<setw(4)<< t_track.GetClusterY_ID(ip);
604     // cout << endl;
605     // cout<< " TRK chi2 "<<setw(10)<<t_track.chi2;
606     // cout <<" | al ";
607     // for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
608     // cout << endl;
609    
610     ///////////////////////////////////////////////////////// FIT START
611     // int ifail=0;
612     // // TrkParams::SetDebugMode();
613     // t_track.FitReset(); //angles set to zero
614     // t_track.Fit(0.,ifail,0,1);//re-evaluate cluster position from level1
615     // if( TMath::IsNaN(t_track.chi2) )continue;
616     // if( ifail != 0 )continue;
617     // t_track.Fit(0.,ifail,0,1);//evaluate again taking into account the track angle
618     // if( TMath::IsNaN(t_track.chi2) )continue;
619     // if( ifail != 0 )continue;
620     ///////////////////////////////////////////////////////// FIT FINISH
621     // cout << "Merged: "<<endl;
622     // cout << "X "; for(int ip=0; ip<6; ip++)cout << t_track.XGood(ip); cout << endl;
623     // cout << "Y "; for(int ip=0; ip<6; ip++)cout << t_track.YGood(ip); cout << endl;
624     //
625     //check if track has been already inserted
626     //
627     multimap<int,TrkTrack>::iterator trkB = trackCandidates[mapIndex].begin();
628     multimap<int,TrkTrack>::iterator trkE = trackCandidates[mapIndex].end();
629     //
630     bool ADD = true;
631     for( multimap<int,TrkTrack>::iterator trkI = trkB; trkI!=trkE; ++trkI){
632     TrkTrack &trackC = (*trkI).second;
633     bool SAME=true;
634     for(int ip=0; ip<6; ip++)
635     if( t_track.GetClusterX_ID(ip) != trackC.GetClusterX_ID(ip) ||
636     t_track.GetClusterY_ID(ip) != trackC.GetClusterY_ID(ip) ||
637     t_track.GetSensor(ip) != trackC.GetSensor(ip) ||
638     false)SAME=false;
639     if(SAME){
640     ADD=false;
641     break;
642     }
643     }
644     // if(_debug)cout <<" ADD ? "<<ADD<<endl;
645     if(ADD){
646    
647     // cout << "-----> INSERT "<<endl;
648     if(_debug){
649     cout << " >> HBM "<<HBM;
650     cout << " X:";for(int ip=0; ip<6; ip++)cout << t_track.XGood(ip);
651     cout << " Y:";for(int ip=0; ip<6; ip++)cout << t_track.YGood(ip);
652     cout << " S:";for(int ip=0; ip<6; ip++)cout << t_track.GetSensor(ip);
653     }
654    
655     ///////////////////////////////////////////////////////// FIT START
656     int ifail=0;
657     t_track.FitReset(); //angles set to zero
658     t_track.Fit(0.,ifail,0,1);//re-evaluate cluster position from level1
659     if( TMath::IsNaN(t_track.chi2) )continue;
660     if( ifail != 0 )continue;
661     t_track.Fit(0.,ifail,0,1);//evaluate again taking into account the track angle
662     if( TMath::IsNaN(t_track.chi2) )continue;
663     if( ifail != 0 )continue;
664     ///////////////////////////////////////////////////////// FIT FINISH
665    
666     if(_debug){
667     cout << " >> al: ";for(int i=0; i<5; i++)cout <<setw(10)<< t_track.al[i];
668     cout << " chi2: "<<setw(10)<< t_track.chi2;
669     cout << endl;
670     }
671    
672     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
673     // evaluated coordinates (to define GF)
674     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
675     int ngf = TrkParams::nGF;
676     float *zgf = TrkParams::zGF;
677     Trajectory tgf = Trajectory(ngf,zgf);
678     tgf.DoTrack(t_track.al);//<<<< integrate the trajectory
679     for(int ip=0; ip<ngf; ip++){
680     t_track.xGF[ip] = tgf.x[ip];
681     t_track.yGF[ip] = tgf.y[ip];
682     }
683     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
684     // evaluated other track info
685     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
686     for(int ip=0; ip<6; ip++){
687     float factor = 0.;
688     factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.axv[ip]/180.),2);
689     factor += TMath::Power(TMath::Tan(TMath::ACos(-1)*t_track.ayv[ip]/180.),2);
690     factor += 1.;
691     factor = TMath::Sqrt(factor);
692     // factor = 1.;///TEMP
693     if(t_track.XGood(ip)){
694     int iclx = t_track.GetClusterX_ID(ip);
695     TrkCluster *clusterx = _trk_l1->GetCluster(iclx);
696     float mip = TrkParams::GetMIP(clusterx->GetLadder()-1,clusterx->view-1);
697     t_track.dedx_x[ip] = clusterx->GetSignal()/(factor>0.?factor:1.)/(mip>0.?mip:1.);
698     t_track.multmaxx[ip] = clusterx->GetMultiplicity()*10000+clusterx->maxs;
699     t_track.seedx[ip] = clusterx->clsignal[clusterx->indmax]/clusterx->clsigma[clusterx->indmax];//clusterx->GetSignalToNoise(0,7);
700     t_track.xpu[ip] = 0.;//complicato....
701     }
702     if(t_track.YGood(ip)){
703     int icly = t_track.GetClusterY_ID(ip);
704     TrkCluster *clustery = _trk_l1->GetCluster(icly);
705     float mip = TrkParams::GetMIP(clustery->GetLadder()-1,clustery->view-1);
706     t_track.dedx_y[ip] = clustery->GetSignal()/(factor>0.?factor:1.)/(mip>0.?mip:1.);
707     t_track.multmaxy[ip] = clustery->GetMultiplicity()*10000+clustery->maxs;
708     t_track.seedy[ip] = clustery->clsignal[clustery->indmax]/clustery->clsigma[clustery->indmax];//clustery->GetSignalToNoise(0,7);
709     t_track.ypu[ip] = 0.;//complicato....
710     }
711     }
712     ///////////////////////////////////////////////////////// end additional info
713    
714     // cout << "ADD "<<endl;
715     trackCandidates[mapIndex].insert(pair<int,TrkTrack>(t_track.GetNX()+t_track.GetNY(), TrkTrack(t_track)));
716     if(t_track.GetNX()+t_track.GetNY()>5 && t_track.chi2<chi2Min) chi2Min=t_track.chi2;
717     iCand++;
718     }//end add-candidate condition
719    
720     }//end first loop on candidates
721     if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
722    
723     //=================================================================
724     // ----------------------------------------------------------------
725     // chi2 selection
726     // ----------------------------------------------------------------
727     //=================================================================
728     pair <multimap<int,TrkTrack>::iterator, multimap<int,TrkTrack>::iterator> nxny[12];
729     for(int inxny=0; inxny<12; inxny++){
730     nxny[inxny] = trackCandidates[mapIndex].equal_range(inxny);
731     for ( multimap<int,TrkTrack>::iterator it = nxny[inxny].first; it!=nxny[inxny].second; ++it){
732     TrkTrack &trackC = (*it).second;
733     // cout << " NX+NY "<< trackC.GetNX()+trackC.GetNY() <<" chi2 "<<trackC.chi2<< " R(GV) "<<trackC.GetRigidity()<<endl;
734     // if( trackC.chi2==0 )continue;
735     if( trackC.GetNX()+trackC.GetNY() > 5 && trackC.chi2>chi2Min )continue;
736     trackCandidates[mapIndex+1].insert(pair<int,TrkTrack>(inxny, TrkTrack(trackC)));
737     }
738     }
739     mapIndex++;
740    
741     if(_debug)cout << "n.candidates "<<trackCandidates[mapIndex].size()<<endl;
742    
743     //=================================================================
744     // ----------------------------------------------------------------
745     // save tracks
746     // ----------------------------------------------------------------
747     //=================================================================
748     trkBegin = trackCandidates[mapIndex].begin();
749     trkEnd = trackCandidates[mapIndex].end();
750     iCand = 0;
751     for( multimap<int,TrkTrack>::iterator trkIter = trkBegin; trkIter!=trkEnd; ++trkIter){
752     TrkTrack trackCand = ((*trkIter).second);
753     trackCand.seqno = iCand;
754     // trackCand.Dump();
755     new(tracks[iCand]) TrkTrack(trackCand);
756     iCand++;
757     };
758    
759    
760    
761     };
762    
763     /**
764     *
765     */
766     void ExtTrkingAlg::Dump(){
767    
768     cout << "|--|--|--|--|--|--|--|--|--|--|--|"<<endl;
769     cout << "ExtTrkingAlg::Dump() "<<endl<<endl;
770     cout << " algorythm "<<setw(10)<<_whichAlg <<endl;
771     cout << " debug "<<setw(10)<< _debug <<endl;
772     cout << " TrkLevel1 "<<setw(10)<<_trk_l1 <<endl;
773     cout << " TrkLevel2 "<<setw(10)<<_trk_l2 <<endl;
774     cout << " CaloLevel1"<<setw(10)<<_cal_l1 <<endl;
775     cout << " CaloLevel2"<<setw(10)<<_cal_l2 <<endl;
776     cout << " ToFLevel2 "<<setw(10)<<_tof_l2 <<endl;
777     cout << "Pre-selection parameters "<<endl;
778     cout << " nClstrMAX "<<setw(10)<<_sel_nClstrMAX <<endl;
779     cout << " nPlaneXMIN "<<setw(10)<<_sel_nPlaneXMIN <<endl;
780     cout << " nPlaneYMIN "<<setw(10)<<_sel_nPlaneYMIN <<endl;
781     cout << "Algorythm parameters "<<endl;
782     cout << " nClFixX "<<setw(10)<<_alg_nClFixX <<endl;
783     cout << " nClFixY "<<setw(10)<<_alg_nClFixY <<endl;
784     cout << " nTrackMAX "<<setw(10)<<_alg_nTrackMAX <<endl;
785     cout << "|--|--|--|--|--|--|--|--|--|--|--|"<<endl;
786    
787     }

  ViewVC Help
Powered by ViewVC 1.1.23