CMSSW Reference Manual

00001 
00002 // File: TrackProducerFP420.cc
00003 // Date: 12.2006
00004 // Description: TrackProducerFP420 for FP420
00005 // Modifications: 
00007 #include "RecoRomanPot/RecoFP420/interface/TrackProducerFP420.h"
00008 #include <stdio.h>
00009 #include <gsl/gsl_fit.h>
00010 #include<vector>
00011 #include <iostream>
00012 using namespace std;
00013 
00014 //#define debugmaxampl
00015 //#define debugvar1
00016 //#define debugvar2
00017 //#define debugvar22
00018 //#define debugsophisticated
00019 //#define debugsophisticated
00020 //#define debug3d
00021 //#define debug3d30
00022 
00023 TrackProducerFP420::TrackProducerFP420(int asn0, int apn0, int azn0, double az420, double azD2, double azD3, double apitchX, double apitchY, double apitchXW, double apitchYW, double aZGapLDet, double aZSiStep, double aZSiPlane, double aZSiDetL, double aZSiDetR, bool aUseHalfPitchShiftInX, bool aUseHalfPitchShiftInY, bool aUseHalfPitchShiftInXW, bool aUseHalfPitchShiftInYW, double adXX, double adYY, float achiCutX, float achiCutY, double azinibeg) {
00024   //
00025   // Everything that depend on the det
00026   //
00027   sn0 = asn0;
00028   pn0 = apn0;
00029   zn0 = azn0;
00030   z420= az420;
00031   zD2 = azD2;
00032   zD3 = azD3;
00033   //zUnit= azUnit;
00034   pitchX = apitchX;
00035   pitchY = apitchY;
00036   pitchXW = apitchXW;
00037   pitchYW = apitchYW;
00038   ZGapLDet = aZGapLDet;
00039   ZSiStep = aZSiStep; 
00040   ZSiPlane = aZSiPlane;
00041   ZSiDetL = aZSiDetL;
00042   ZSiDetR = aZSiDetR;
00043   UseHalfPitchShiftInX = aUseHalfPitchShiftInX;
00044   UseHalfPitchShiftInY = aUseHalfPitchShiftInY;
00045   UseHalfPitchShiftInXW = aUseHalfPitchShiftInXW;
00046   UseHalfPitchShiftInYW = aUseHalfPitchShiftInYW;
00047   dXX = adXX;
00048   dYY = adYY;
00049   chiCutX = achiCutX;
00050   chiCutY = achiCutY;
00051   zinibeg = azinibeg;
00052 
00053 #ifdef debugmaxampl
00054   std::cout << "TrackProducerFP420: call constructor" << std::endl;
00055   std::cout << " sn0= " << sn0 << " pn0= " << pn0 << " zn0= " << zn0 << std::endl;
00056   std::cout << " zD2= " << zD2 << " zD3= " << zD3 << " zinibeg= " << zinibeg << std::endl;
00057   //std::cout << " zUnit= " << zUnit << std::endl;
00058   std::cout << " pitchX= " << pitchX << " pitchY= " << pitchY << std::endl;
00059   std::cout << " ZGapLDet= " << ZGapLDet << std::endl;
00060   std::cout << " ZSiStep= " << ZSiStep << " ZSiPlane= " << ZSiPlane << std::endl;
00061   std::cout << " ZSiDetL= " <<ZSiDetL  << " ZSiDetR= " << ZSiDetR << std::endl;
00062   std::cout << " UseHalfPitchShiftInX= " << UseHalfPitchShiftInX << " UseHalfPitchShiftInY= " << UseHalfPitchShiftInY << std::endl;
00063   std::cout << "TrackProducerFP420:----------------------" << std::endl;
00064   std::cout << " dXX= " << dXX << " dYY= " << dYY << std::endl;
00065   std::cout << " chiCutX= " << chiCutX << " chiCutY= " << chiCutY << std::endl;
00066 #endif
00067 }
00069 
00070 
00071 
00072 
00073 
00075 std::vector<TrackFP420> TrackProducerFP420::trackFinderSophisticated(edm::Handle<ClusterCollectionFP420> input, int det){
00077   
00078   std::vector<TrackFP420> rhits;
00079   int restracks = 10;// max # tracks
00080   rhits.reserve(restracks); 
00081   rhits.clear();
00082   double Ax[10]; double Bx[10]; double Cx[10]; int Mx[10];
00083   double Ay[10]; double By[10]; double Cy[10]; int My[10];
00084   double AxW[10]; double BxW[10]; double CxW[10]; int MxW[10];
00085   double AyW[10]; double ByW[10]; double CyW[10]; int MyW[10];
00086 #ifdef debugsophisticated
00087   std::cout << "TrackProducerFP420: Start trackFinderSophisticated " << std::endl; 
00088 #endif
00091 // zn0 is the same as xytype
00092   if( zn0 < 1 || zn0 > 4 ){
00093     std::cout << "TrackProducerFP420:ERROR in trackFinderSophisticated: check zn0 (xytype) = " << zn0 << std::endl; 
00094     return rhits;
00095   }
00096 // sn0= 3 - 2St configuration, sn0= 4 - 3St configuration 
00097   if( sn0 < 3 || zn0 > 4 ){
00098     std::cout << "TrackProducerFP420:ERROR in trackFinderSophisticated: check sn0 (configuration) = " << sn0 << std::endl; 
00099     return rhits;
00100   }
00102   int zbeg = 1, zmax=3;// XY
00103   //  if( zn0==1){
00104   //              zmax=2; // Y
00105   //  }
00106   //  else if( zn0==2){
00107   //    zbeg = 2, zmax=3; // X
00108   // }
00110   //   .
00111   int reshits = 15;// max # cl for every X and Y plane
00112   //  int resplanes = 30;
00113   int nX[30], nY[30];// resplanes =30 NUMBER OF PLANES; nX, nY - # cl for every X and Y plane
00114   int uX[30], uY[30];// resplanes =30 NUMBER OF PLANES; nX, nY - current # cl used for every X and Y plane
00115   double zX[15][30], xX[15][30], wX[15][30];
00116   double zY[15][30], yY[15][30], wY[15][30];
00117   double             yXW[15][30], wXW[15][30];
00118   double             xYW[15][30], wYW[15][30];
00119   bool qX[15][30], qY[15][30];
00120   //   .
00121   int txf = 0; int txs1 = 0; int txss = 0;
00122   int tyf = 0; int tys1 = 0; int tyss = 0;
00123   //   .
00124   double pitch=0.;
00125   double pitchW=0.;
00126   if(zn0==1){
00127     pitch=pitchY;
00128     pitchW=pitchYW;
00129   }
00130   else if(zn0==2){
00131     pitch=pitchX;
00132     pitchW=pitchXW;
00133   }
00134 
00135 
00136      //current change of geometry:
00137     float Xshift = pitch/2.;
00138     float Yshift = pitchW/2.;
00139     
00140     //
00141     
00142     for (int sector=1; sector < sn0; sector++) {
00143       for (int zmodule=1; zmodule<pn0; zmodule++) {
00144         for (int zside=zbeg; zside<zmax; zside++) {
00145           
00146           // index iu is a continues numbering of 3D detector of FP420 (detector ID)
00147           int sScale = 2*(pn0-1), dScale = 2*(pn0-1)*(sn0-1);
00148           int zScale=2;  unsigned int iu = dScale*(det - 1)+sScale*(sector - 1)+zScale*(zmodule - 1)+zside;
00149           //    unsigned int ii = sScale*(sector - 1)/2 + (zmodule - 1) + 1;
00150 
00151           //      unsigned int ii = sScale*(sector - 1)/2 + (zmodule - 1) ; // 0-19   --> 20 items
00152           unsigned int ii = iu-1-dScale*(det - 1);// 0-29   --> 30 items
00153           
00154           double kplane = -(pn0-1)/2 - 0.5  +  (zmodule-1); 
00155           
00156           
00157           double zdiststat = 0.;
00158           if(sn0<4) {
00159             if(sector==2) zdiststat = zD3;
00160           }
00161           else {
00162             if(sector==2) zdiststat = zD2;
00163             if(sector==3) zdiststat = zD3;
00164           }
00165           double zcurrent = zinibeg + z420 + (ZSiStep-ZSiPlane)/2  + kplane*ZSiStep + zdiststat;  
00166           //double zcurrent = zinibeg +(ZSiStep-ZSiPlane)/2  + kplane*ZSiStep + (sector-1)*zUnit;  
00167           
00168           if(zside==1){
00169             zcurrent += (ZGapLDet+ZSiDetL/2);
00170           }
00171           if(zside==2){
00172             zcurrent += (ZGapLDet+ZSiDetR/2)+ZSiPlane/2;
00173           }
00174           //   .
00175           //
00176           if(det == 2) zcurrent = -zcurrent;
00177           //
00178           //
00179           //   .
00180           // local - global systems with possible shift of every second plate:
00181           
00182           // for zn0=1
00183           float dYYcur = dYY;// XSiDet/2.
00184           float dYYWcur = dXX;//(BoxYshft+dYGap) + (YSi - YSiDet)/2. = 12.7
00185           // for zn0=2
00186           float dXXcur = dXX;//(BoxYshft+dYGap) + (YSi - YSiDet)/2. = 12.7
00187           float dXXWcur = dYY;// XSiDet/2.
00188           //   .
00189           if(zside==2) {
00190             // X-type: x-coord
00191             if (UseHalfPitchShiftInX == true){
00192               dXXcur += Xshift;
00193             }
00194             // X-type: y-coord
00195             if (UseHalfPitchShiftInXW == true){
00196               dXXWcur -= Yshift;
00197             }
00198           }
00199           //
00200           
00201           //   .
00202           //   GET CLUSTER collection  !!!!
00203           //   .
00204           //============================================================================================================ put into currentclust
00205           std::vector<ClusterFP420> currentclust;
00206           currentclust.clear();
00207           ClusterCollectionFP420::Range outputRange;
00208           outputRange = input->get(iu);
00209           // fill output in currentclust vector (for may be sorting? or other checks)
00210           ClusterCollectionFP420::ContainerIterator sort_begin = outputRange.first;
00211           ClusterCollectionFP420::ContainerIterator sort_end = outputRange.second;
00212           for ( ;sort_begin != sort_end; ++sort_begin ) {
00213             //  std::sort(currentclust.begin(),currentclust.end());
00214             currentclust.push_back(*sort_begin);
00215           } // for
00216           
00217 #ifdef debugsophisticated
00218           std::cout << "TrackProducerFP420: currentclust.size = " << currentclust.size() << std::endl; 
00219 #endif
00220           //============================================================================================================
00221           
00222           vector<ClusterFP420>::const_iterator simHitIter = currentclust.begin();
00223           vector<ClusterFP420>::const_iterator simHitIterEnd = currentclust.end();
00224           
00225           if(zn0 ==1){
00226             nY[ii] = 0;// # cl in every Y plane (max is reshits)
00227             uY[ii] = 0;// current used # cl in every X plane 
00228           }
00229           else if(zn0 ==2){
00230             nX[ii] = 0;// # cl in every X plane (max is reshits)
00231             uX[ii] = 0;// current used # cl in every X plane 
00232           }
00233           // loop in #clusters
00234           for (;simHitIter != simHitIterEnd; ++simHitIter) {
00235             const ClusterFP420 icluster = *simHitIter;
00236             
00237             // fill vectors for track reconstruction
00238             
00239             
00240             //disentangle complicated pattern recognition of hits?
00241             // Y:
00242             if(zn0 ==1){
00243               nY[ii]++;         
00244               if(nY[ii]>reshits){
00245                 nY[ii]=reshits;
00246                 std::cout << "WARNING-ERROR:TrackproducerFP420: currentclust.size()= " << currentclust.size() <<" bigger reservated number of hits" << " zcurrent=" << zY[nY[ii]-1][ii] << " ii= "  << ii << std::endl;
00247               }
00248               zY[nY[ii]-1][ii] = zcurrent;
00249               yY[nY[ii]-1][ii] = icluster.barycenter()*pitch;
00250               xYW[nY[ii]-1][ii] = icluster.barycenterW()*pitchW;
00251               // go to global system:
00252               yY[nY[ii]-1][ii] = yY[nY[ii]-1][ii] - dYYcur; 
00253               wY[nY[ii]-1][ii] = 1./(icluster.barycerror()*pitch);//reciprocal of the variance for each datapoint in y
00254               wY[nY[ii]-1][ii] *= wY[nY[ii]-1][ii];//reciprocal of the variance for each datapoint in y
00255               xYW[nY[ii]-1][ii] =-(xYW[nY[ii]-1][ii]+dYYWcur); 
00256               wYW[nY[ii]-1][ii] = 1./(icluster.barycerrorW()*pitchW);//reciprocal of the variance for each datapoint in y
00257               wYW[nY[ii]-1][ii] *= wYW[nY[ii]-1][ii];//reciprocal of the variance for each datapoint in y
00258               qY[nY[ii]-1][ii] = true;
00259               if(nY[ii]==reshits) break;
00260             }
00261             // X:
00262             else if(zn0 ==2){
00263               nX[ii]++; 
00264               if(nX[ii]>reshits){
00265                 std::cout << "WARNING-ERROR:TrackproducerFP420: currentclust.size()= " << currentclust.size() <<" bigger reservated number of hits" << std::endl;
00266                 nX[ii]=reshits;
00267               }
00268               zX[nX[ii]-1][ii] = zcurrent;
00269               xX[nX[ii]-1][ii] = icluster.barycenter()*pitch;
00270               yXW[nX[ii]-1][ii] = icluster.barycenterW()*pitchW;
00271               // go to global system:
00272               xX[nX[ii]-1][ii] =-(xX[nX[ii]-1][ii]+dXXcur); 
00273               wX[nX[ii]-1][ii] = 1./(icluster.barycerror()*pitch);//reciprocal of the variance for each datapoint in y
00274               wX[nX[ii]-1][ii] *= wX[nX[ii]-1][ii];//reciprocal of the variance for each datapoint in y
00275               yXW[nX[ii]-1][ii] = yXW[nX[ii]-1][ii] - dXXWcur; 
00276               wXW[nX[ii]-1][ii] = 1./(icluster.barycerrorW()*pitchW);//reciprocal of the variance for each datapoint in y
00277               wXW[nX[ii]-1][ii] *= wXW[nX[ii]-1][ii];//reciprocal of the variance for each datapoint in y
00278               qX[nX[ii]-1][ii] = true;
00279 #ifdef debugsophisticated
00280               std::cout << "trackFinderSophisticated: nX[ii]= " << nX[ii]<< " ii = " << ii << " zcurrent = " << zcurrent << " xX[nX[ii]-1][ii] = " << xX[nX[ii]-1][ii] << std::endl;
00281               std::cout << " wX[nX[ii]-1][ii] = " << wX[nX[ii]-1][ii] << " wXW[nX[ii]-1][ii] = " << wXW[nX[ii]-1][ii] << std::endl;
00282               std::cout << " -icluster.barycenter()*pitch = " << -icluster.barycenter()*pitch << " -dXXcur = " << -dXXcur << std::endl;
00283               std::cout << "============================================================" << std::endl;
00284 #endif
00285               if(nX[ii]==reshits) break;
00286             }
00287             
00288           } // for loop in #clusters (can be breaked)
00289           
00290           // Y:
00291           if(zn0 ==1){
00292             if(nY[ii] != 0) {  /* # Y-planes w/ clusters */
00293               ++tyf; if(sector==1) ++tys1; if(sector==(sn0-1)) ++tyss;
00294             }     
00295           }
00296           // X:
00297           else if(zn0 ==2){
00298             if(nX[ii] != 0) {  /* # X-planes w/ clusters */
00299               ++txf; if(sector==1) ++txs1; if(sector==(sn0-1)) ++txss;
00300             }     
00301           }
00302           //================================== end of for loops in continuius number iu:
00303         }   // for zside
00304       }   // for zmodule
00305     }   // for sector
00306 #ifdef debugsophisticated
00307     std::cout << "trackFinderSophisticated: tyf= " << tyf<< " tys1 = " << tys1 << " tyss = " << tyss << std::endl;
00308     std::cout << "trackFinderSophisticated: txf= " << txf<< " txs1 = " << txs1 << " txss = " << txss << std::endl;
00309     std::cout << "============================================================" << std::endl;
00310 #endif
00311     
00312     //===========================================================================================================================
00313     //===========================================================================================================================
00314     //===========================================================================================================================
00315     //======================    start road finder   =============================================================================
00316     //===========================================================================================================================
00317 
00318   //  int nitMax=5;// max # iterations to find track
00319   int nitMax=5;// max # iterations to find track
00320 
00321   // criteria for track selection: 
00322   // track is selected if for 1st station #cl >=pys1Cut
00323   //  int  pys1Cut = 5, pyssCut = 5, pyallCut=12;
00324   //  int  pys1Cut = 1, pyssCut = 1, pyallCut= 3;
00325   int  pys1Cut = 3, pyssCut = 3, pyallCut= 6;
00326   
00327   //  double yyyvtx = 0.0, xxxvtx = -15;  //mm
00328   
00330   //
00331   // for configuration: 3St, 1m for 1-2 St:
00332   // double sigman=0.1, ssigma = 1.0, sigmam=0.15;/* ssigma is foreseen to match 1st point of 2nd Station*/
00333   //
00334   // for equidistant 3 Stations:
00335   //
00336   // for tests:
00337   //  double sigman=118., ssigma = 299., sigmam=118.;
00338   // RMS1=0.013, RMS2 = 1.0, RMS3 = 0.018 see plots d1XCL, d2XCL, d3XCL
00339   //
00340   //  double sigman=0.05, ssigma = 2.5, sigmam=0.06;
00341   //  double sigman=0.18, ssigma = 1.8, sigmam=0.18;
00342   //  double sigman=0.18, ssigma = 2.9, sigmam=0.18;
00343   //
00345   // for 3 Stations:
00346   // LAST:
00347   double sigman=0.18, ssigma = 2.5, sigmam=0.18;
00348   if( sn0 < 4 ){
00349     // for 2 Stations:
00350     // sigman=0.24, ssigma = 4.2, sigmam=0.33;
00351     //  sigman=0.18, ssigma = 3.9, sigmam=0.18;
00352     // sigman=0.18, ssigma = 3.6, sigmam=0.18;
00353     sigman=0.18, ssigma = 3.3, sigmam=0.18;
00354   }
00355 #ifdef debugsophisticated
00356   std::cout << "trackFinderSophisticated: ssigma= " << ssigma << std::endl;
00357 #endif
00358 
00360   
00361   /* ssigma = 3. * 8000.*(0.025+0.009)/sqrt(pn0-1)/100. = 2.9 mm(!!!)----
00362      ssigma is reduced by factor k_reduced = (sn0-1)-sector+1 = sn0-sector
00363      # Stations  currentStation
00364      2Stations:     sector=2,         sn0=3 , sn0-sector = 1 --> k_reduced = 1
00365      3Stations:     sector=2,         sn0=4 , sn0-sector = 2 --> k_reduced = 2
00366      3Stations:     sector=3,         sn0=4 , sn0-sector = 1 --> k_reduced = 1
00367   */
00368   int numberXtracks=0, numberYtracks=0, totpl = 2*(pn0-1)*(sn0-1); double sigma;
00369 
00370   //  for (int zside=zbeg; zside<zmax; ++zside) {
00371   for (int zsidezn0=zn0; zsidezn0<zn0+1; ++zsidezn0) {
00372 #ifdef debugsophisticated
00373   std::cout << "trackFinderSophisticated: zsidezn0= " << zsidezn0 << std::endl;
00374 #endif
00375 
00376     //
00377     //
00378     double tg0 = 0.;
00379     int qAcl[30], qAii[30], fip=0, niteration = 0;
00380     int ry = 0, rys1 = 0, ryss = 0;
00381     int tas1=tys1, tass=tyss, taf=tyf;
00382     bool SelectTracks = true;
00383     //
00385   //   .
00386 
00387   double yA[15][30], zA[15][30], wA[15][30]; int nA[30], uA[30]; bool qA[15][30];
00388     //
00389     // Y:
00390   //======================    start road finder  for zsidezn0 = 1      ===========================================================
00391     if(zsidezn0 ==1){
00392   //===========================================================================================================================
00393       numberYtracks=0;  
00394       tg0= 3*1./(800.+20.); // for Y: 1cm/...   *3 - 3sigma range
00395       tas1=tys1;
00396       tass=tyss;
00397       taf=tyf;
00398       for (int ii=0; ii < totpl; ++ii) {
00399 #ifdef debugsophisticated
00400   std::cout << "trackFinderSophisticated: ii= " << ii << " nY[ii]= " << nY[ii] << std::endl;
00401   std::cout << "trackFinderSophisticated: ii= " << ii << " uY[ii]= " << uY[ii] << std::endl;
00402 #endif
00403         nA[ii] = nY[ii];
00404         uA[ii] = uY[ii];
00405         for (int cl=0; cl<nA[ii]; ++cl) {
00406 #ifdef debugsophisticated
00407   std::cout << " cl= " << cl << " yY[cl][ii]= " << yY[cl][ii] << std::endl;
00408   std::cout << " zY[cl][ii]= " << zY[cl][ii] << " wY[cl][ii]= " << wY[cl][ii] << " qY[cl][ii]= " << qY[cl][ii] << std::endl;
00409 #endif
00410           yA[cl][ii] = yY[cl][ii];
00411           zA[cl][ii] = zY[cl][ii];
00412           wA[cl][ii] = wY[cl][ii];
00413           qA[cl][ii] = qY[cl][ii];
00414         }
00415       }
00416   //===========================================================================================================================
00417     }// if zsidezn0 ==1
00418     // X:
00419   //======================    start road finder  for zside = 2      ===========================================================
00420     else if(zsidezn0 ==2){
00421   //===========================================================================================================================
00422       numberXtracks=0;  
00423       tg0= 3*2./(800.+20.); // for X: 2cm/...   *3 - 3sigma range
00424       tas1=txs1;
00425       tass=txss;
00426       taf=txf;
00427       for (int ii=0; ii < totpl; ++ii) {
00428 #ifdef debugsophisticated
00429   std::cout << "trackFinderSophisticated: ii= " << ii << " nX[ii]= " << nX[ii] << std::endl;
00430   std::cout << "trackFinderSophisticated: ii= " << ii << " uX[ii]= " << uX[ii] << std::endl;
00431 #endif
00432         nA[ii] = nX[ii];
00433         uA[ii] = uX[ii];
00434         for (int cl=0; cl<nA[ii]; ++cl) {
00435 #ifdef debugsophisticated
00436   std::cout << " cl= " << cl << " xX[cl][ii]= " << xX[cl][ii] << std::endl;
00437   std::cout << " zX[cl][ii]= " << zX[cl][ii] << " wX[cl][ii]= " << wX[cl][ii] << " qX[cl][ii]= " << qX[cl][ii] << std::endl;
00438 #endif
00439           yA[cl][ii] = xX[cl][ii];
00440           zA[cl][ii] = zX[cl][ii];
00441           wA[cl][ii] = wX[cl][ii];
00442           qA[cl][ii] = qX[cl][ii];
00443         }
00444       }
00445   //===========================================================================================================================
00446     }// if zsidezn0 ==zn0
00447 
00448 
00449     
00450   //======================    start road finder        ====================================================
00451     do {
00452       double fyY[30], fzY[30], fwY[30];
00453       double fyYW[30],         fwYW[30];
00454       int py = 0, pys1 = 0, pyss = 0;
00455       bool NewStation = false, py1first = false;
00456       for (int sector=1; sector < sn0; ++sector) {
00457         double tav=0., t1=0., t2=0., t=0., sm;
00458         int stattimes=0;
00459         if( sector != 1 ) {
00460           NewStation = true;  
00461         }
00462         for (int zmodule=1; zmodule<pn0; ++zmodule) {
00463           for (int zside=zbeg; zside<zmax; zside++) {
00464             
00465             // index iu is a continues numbering of 3D detector of FP420 (detector ID)
00466             int sScale = 2*(pn0-1);
00467             int zScale=2;  unsigned int iu = sScale*(sector - 1)+zScale*(zmodule - 1)+zside;
00468             unsigned int ii = iu-1;// 0-29   --> 30 items
00469             
00470             
00471             if(nA[ii]!=0  && uA[ii]!= nA[ii]) { 
00472               
00473               ++py; if(sector==1) ++pys1; if(sector==(sn0-1)) ++pyss;
00474               if(py==2 && sector==1) { 
00475                 double dymin=9999999., df2; int cl2=-1;
00476                 for (int cl=0; cl<nA[ii]; ++cl) {
00477                   if(qA[cl][ii]){
00478                     df2 = abs(fyY[fip]-yA[cl][ii]);
00479                     if(df2 < dymin) {
00480                       dymin = df2;
00481                       cl2=cl;
00482                     }//if(df2           
00483                   }//if(qA              
00484                 }//for(cl
00485                 if(cl2!=-1){
00486                   t=(yA[cl2][ii]-fyY[fip])/(zA[cl2][ii]-fzY[fip]);
00487                   t1 = t*wA[cl2][ii];
00488                   t2 = wA[cl2][ii];
00489 #ifdef debugsophisticated
00490                   std::cout << " t= " << t << " tg0= " << tg0 << std::endl;
00491 #endif
00492                   if(abs(t)<tg0) { 
00493                     qA[cl2][ii] = false;//point is taken, mark it for not using again
00494                     fyY[py-1]=yA[cl2][ii];
00495                     fzY[py-1]=zA[cl2][ii];
00496                     fwY[py-1]=wA[cl2][ii];
00497                     qAcl[py-1] = cl2;
00498                     qAii[py-1] = ii;
00499                     ++uA[ii];
00500 #ifdef debugsophisticated
00501                     std::cout << " point is taken, mark it for not using again uA[ii]= " << uA[ii] << std::endl;
00502 #endif
00503                     if(uA[ii]==nA[ii]){/* no points anymore for this plane */
00504                       ++ry; if(sector==1) ++rys1; if(sector==(sn0-1)) ++ryss;
00505                     }//if(uA
00506                   }//if abs
00507                   else{
00508                     py--; if(sector==1) pys1--;  if(sector==(sn0-1)) pyss--;
00509                     t1 -= t*wA[cl2][ii]; t2 -= wA[cl2][ii];
00510                   }//if(abs
00511                 }//if(cl2!=-1
00512                 else{
00513                   py--; if(sector==1) pys1--;  if(sector==(sn0-1)) pyss--;
00514                 }//if(cl2!=-1
00515               }//if(py==2
00516               else {
00517                 int clcurr=-1;
00518                 for (int cl=0; cl<nA[ii]; ++cl) {
00519                   if(qA[cl][ii]){
00520                     clcurr = cl;
00521                     if(py<3 ){
00522                       if(py==1) { 
00523                         py1first = true;
00524                         fip=py-1;
00525                         qA[cl][ii] = false;//point is taken, mark it for not using again
00526                         fyY[py-1]=yA[cl][ii];
00527                         fzY[py-1]=zA[cl][ii];
00528                         fwY[py-1]=wA[cl][ii];
00529                         qAcl[py-1] = cl;
00530                         qAii[py-1] = ii;
00531                         ++uA[ii];
00532 #ifdef debugsophisticated
00533                         std::cout << " point is taken, mark it uA[ii]= " << uA[ii] << std::endl;
00534 #endif
00535                       }//if py=1
00536                       if(uA[ii]==nA[ii]){/* no points anymore for this plane */
00537                         ++ry; if(sector==1) ++rys1; if(sector==(sn0-1)) ++ryss;
00538                       }//if(uA
00539                     }//py<3
00540                     else {
00541                       if(NewStation){
00542                         if( sn0 < 4 ) {
00543                           sigma = ssigma;
00544                         }
00545                         else {
00546                           sigma = ssigma/(sn0-1-sector);
00547                         }
00548                         //      std::cout << " sector= " << sector << " sn0= " << sn0 << " sigma= " << sigma << std::endl;
00549                         //      std::cout << " stattimes= " << stattimes << " ssigma= " << ssigma << " sigmam= " << sigmam << std::endl;
00550 
00551                         //sigma = ssigma/(sn0-sector);
00552                         //if(stattimes==1 || sector==3 ) sigma = msigma * sqrt(1./wA[cl][ii]);
00553 
00554                         if(stattimes==1 || sector==3 ) sigma = sigmam; // (1st $3rd Stations for 3St. configur. ), 1st only for 2St. conf.
00555                         //      if(stattimes==1 || sector==(sn0-1) ) sigma = sigmam;
00556 
00557                         double cov00, cov01, cov11, c0Y, c1Y, chisqY;
00558                         gsl_fit_wlinear (fzY, 1, fwY, 1, fyY, 1, py-1, 
00559                                          &c0Y, &c1Y, &cov00, &cov01, &cov11, 
00560                                          &chisqY);
00561                         sm = c0Y+ c1Y*zA[cl][ii];
00562                         
00563 #ifdef debugsophisticated
00564                           std::cout << " sector= " << sector << " sn0= " << sn0 << " sigma= " << sigma << std::endl;
00565                           std::cout << " stattimes= " << stattimes << " ssigma= " << ssigma << " sigmam= " << sigmam << std::endl;
00566                           std::cout << " sm= " << sm << " c0Y= " << c0Y << " c1Y= " << c1Y << " chisqY= " << chisqY << std::endl;
00567                           std::cout << " zA[cl][ii]= " << zA[cl][ii] << " ii= " << ii << " cl= " << cl << std::endl;
00568                         for (int ct=0; ct<py-1; ++ct) {
00569                           std::cout << " py-1= " << py-1 << " fzY[ct]= " << fzY[ct] << std::endl;
00570                           std::cout << " fyY[ct]= " << fyY[ct] << " fwY[ct]= " << fwY[ct] << std::endl;
00571                         }
00572 #endif
00573                         
00574                       }//NewStation 1
00575                       else{
00576                         t=(yA[cl][ii]-fyY[fip])/(zA[cl][ii]-fzY[fip]);
00577                         t1 += t*wA[cl][ii];
00578                         t2 += wA[cl][ii];
00579                         tav=t1/t2;
00580                         sm = fyY[fip]+tav*(zA[cl][ii]-fzY[fip]);
00581                         //sigma = nsigma * sqrt(1./wA[cl][ii]);
00582                         sigma = sigman;
00583                       }
00584 
00585                       double diffpo = yA[cl][ii]-sm;
00586 #ifdef debugsophisticated
00587                           std::cout << " diffpo= " << diffpo << " yA[cl][ii]= " << yA[cl][ii] << " sm= " << sm << " sigma= " << sigma << std::endl;
00588 #endif
00589                       
00590                       if(abs(diffpo) < sigma ) {
00591                         if(NewStation){
00592                           ++stattimes;
00593                           if(stattimes==1) {
00594                             fip=py-1;
00595                             t1 = 0; t2 = 0;
00596                           }
00597                           else if(stattimes==2){
00598                             NewStation = false; 
00599                             t=(yA[cl][ii]-fyY[fip])/(zA[cl][ii]-fzY[fip]);
00600                             //t1 += t*wA[cl][ii];
00601                             //t2 += wA[cl][ii];
00602                             t1 = t*wA[cl][ii];
00603                             t2 = wA[cl][ii];
00604                           }//if(stattime
00605                         }//if(NewStation 2
00606                         fyY[py-1]=yA[cl][ii];
00607                         fzY[py-1]=zA[cl][ii];
00608                         fwY[py-1]=wA[cl][ii];
00609                         qA[cl][ii] = false;//point is taken, mark it for not using again
00610                         qAcl[py-1] = cl;
00611                         qAii[py-1] = ii;
00612                         ++uA[ii];
00613 #ifdef debugsophisticated
00614                         std::cout << " 3333 point is taken, mark it uA[ii]= " << uA[ii] << std::endl;
00615 #endif
00616                         if(uA[ii]==nA[ii]){/* no points anymore for this plane */
00617                           ++ry; if(sector==1) ++rys1; if(sector==(sn0-1)) ++ryss;
00618                         }//if(cl==
00619                         //  break; // to go on neyt plane
00620                       }//if abs
00621                       else{
00622                         t1 -= t*wA[cl][ii]; t2 -= wA[cl][ii];
00623                       }//if abs
00624                     }// if py<3 and else py>3
00625                     
00626                     if(!qA[cl][ii]) break;// go on neyt plane if point is found among clusters of current plane;
00627                   }// if qA
00628                 }// for cl     --  can be break and return to "for zmodule"
00629                 if( (py!=1 && clcurr != -1 && qA[clcurr][ii]) || (py==1 && !py1first)) { 
00630                   // if point is not found - continue natural loop, but reduce py 
00631                   py--; if(sector==1) pys1--;  if(sector==(sn0-1)) pyss--;
00632                 }//if(py!=1
00633               }//if(py==2 else 
00634             }//if(nA !=0           : inside  this if( -  ask  ++py
00635           }// for zside
00636         }// for zmodule
00637       }// for sector
00638       //============
00639       
00640       
00641 #ifdef debugsophisticated
00642       std::cout << "END: pys1= " << pys1 << " pyss = " << pyss << " py = " << py << std::endl;
00643 #endif
00644       // apply criteria for track selection: 
00645       // do not take track if 
00646       if( pys1 < pys1Cut || pyss < pyssCut || py < pyallCut ){
00647         //      if( pys1<3 || pyss<2 || py<4 ){
00648       }
00649       // do fit:
00650       else{
00652         double cov00, cov01, cov11;
00653         double c0Y, c1Y, chisqY;
00654         gsl_fit_wlinear (fzY, 1, fwY, 1, fyY, 1, py, 
00655                          &c0Y, &c1Y, &cov00, &cov01, &cov11, 
00656                          &chisqY);
00658 #ifdef debugsophisticated
00659         float chindfx;
00660         if(py>2) {
00661           chindfx = chisqY/(py-2);
00662         }
00663         else{
00664           //      chindfy = chisqY;
00665           chindfx = 9999;
00666         }//py
00667         std::cout << " Do FIT XZ: chindfx= " << chindfx << std::endl;
00668 #endif
00669 
00671 #ifdef debugsophisticated
00672         std::cout << " preparation for second order fit for Wide pixels= " << std::endl;
00673 #endif
00674         for (int ipy=0; ipy<py; ++ipy) {
00675           if(zsidezn0 ==1){
00676             fyYW[ipy]=xYW[qAcl[ipy]][qAii[ipy]];
00677             fwYW[ipy]=wYW[qAcl[ipy]][qAii[ipy]];
00678 #ifdef debugsophisticated
00679         std::cout << " ipy= " << ipy << std::endl;
00680         std::cout << " qAcl[ipy]= " << qAcl[ipy] << " qAii[ipy]= " << qAii[ipy] << std::endl;
00681         std::cout << " fyYW[ipy]= " << fyYW[ipy] << " fwYW[ipy]= " << fwYW[ipy] << std::endl;
00682 #endif
00683           }
00684           else if(zsidezn0 ==2){
00685             fyYW[ipy]=yXW[qAcl[ipy]][qAii[ipy]];
00686             fwYW[ipy]=wXW[qAcl[ipy]][qAii[ipy]];
00687 #ifdef debugsophisticated
00688         std::cout << " ipy= " << ipy << std::endl;
00689         std::cout << " qAcl[ipy]= " << qAcl[ipy] << " qAii[ipy]= " << qAii[ipy] << std::endl;
00690         std::cout << " fyYW[ipy]= " << fyYW[ipy] << " fwYW[ipy]= " << fwYW[ipy] << std::endl;
00691 #endif
00692           }
00693         }
00694 #ifdef debugsophisticated
00695         std::cout << " start second order fit for Wide pixels= " << std::endl;
00696 #endif
00697         double wov00, wov01, wov11;
00698         double w0Y, w1Y, whisqY;
00699         gsl_fit_wlinear (fzY, 1, fwYW, 1, fyYW, 1, py, 
00700                          &w0Y, &w1Y, &wov00, &wov01, &wov11, 
00701                          &whisqY);
00703         float chindfy;
00704         if(py>2) {
00705           chindfy = chisqY/(py-2);
00706         }
00707         else{
00708           //      chindfy = chisqY;
00709           chindfy = 9999;
00710         }//py
00711         
00712 #ifdef debugsophisticated
00713         std::cout << " chindfy= " << chindfy << " chiCutY= " << chiCutY << std::endl;
00714 #endif
00715         if(zsidezn0 ==1){
00716           if(chindfy < chiCutX ) {
00717             ++numberYtracks;
00718             Ay[numberYtracks-1] = c0Y; 
00719             By[numberYtracks-1] = c1Y; 
00720             Cy[numberYtracks-1] = chisqY; 
00721             My[numberYtracks-1] = py;
00722             AyW[numberYtracks-1] = w0Y; 
00723             ByW[numberYtracks-1] = w1Y; 
00724             CyW[numberYtracks-1] = whisqY; 
00725             MyW[numberYtracks-1] = py;
00726 #ifdef debugsophisticated
00727             if(py>30) {
00728               std::cout << " niteration = " << niteration << std::endl;
00729               std::cout << " chindfy= " << chindfy << " py= " << py << std::endl;
00730               std::cout << " c0Y= " << c0Y << " c1Y= " << c1Y << std::endl;
00731               std::cout << " pys1= " << pys1 << " pyss = " << pyss << std::endl;
00732             }
00733 #endif
00734           }//chindfy
00735         }
00736         else if(zsidezn0 ==2){
00737           if(chindfy < chiCutY ) {
00738             ++numberXtracks;
00739             Ax[numberXtracks-1] = c0Y; 
00740             Bx[numberXtracks-1] = c1Y; 
00741             Cx[numberXtracks-1] = chisqY; 
00742             Mx[numberXtracks-1] = py;
00743             AxW[numberXtracks-1] = w0Y; 
00744             BxW[numberXtracks-1] = w1Y; 
00745             CxW[numberXtracks-1] = whisqY; 
00746             MxW[numberXtracks-1] = py;
00747 #ifdef debugsophisticated
00748               std::cout << " niteration = " << niteration << std::endl;
00749               std::cout << " chindfx= " << chindfy << " px= " << py << std::endl;
00750               std::cout << " c0X= " << c0Y << " c1X= " << c1Y << std::endl;
00751               std::cout << " pxs1= " << pys1 << " pxss = " << pyss << std::endl;
00752 #endif
00753           }//chindfy
00754         }
00755         
00756         
00757       }//  if else
00758         
00759       // do not select tracks anymore if
00760 #ifdef debugsophisticated
00761       std::cout << " numberYtracks= " << numberYtracks << std::endl;
00762       std::cout << " numberXtracks= " << numberXtracks << std::endl;
00763       std::cout << " pys1= " << pys1 << " pyss = " << pyss << " py = " << py << std::endl;
00764       std::cout << " tas1= " << tas1 << " tass = " << tass << " taf = " << taf << std::endl;
00765       std::cout << " rys1= " << rys1 << " ryss = " << ryss << " ry = " << ry << std::endl;
00766       std::cout << " tas1-rys1= " << tas1-rys1 << " tass-ryss = " << tass-ryss << " taf-ry = " << taf-ry << std::endl;
00767       std::cout << "---------------------------------------------------------- " << std::endl;
00768 #endif
00769       // let's decide: do we continue track finder procedure
00770       if( tas1-rys1<pys1Cut || tass-ryss<pyssCut || taf-ry<pyallCut  ){
00771         SelectTracks = false;
00772       }
00773       else{
00774         ++niteration;
00775 #ifdef debugsophisticated
00776         if(niteration > nitMax-1){
00777           std::cout << "Neyt iteration, niteration >= " << niteration << std::endl;
00778         }
00779 #endif
00780       }
00781       
00782     } while(SelectTracks && niteration < nitMax );      
00783   //======================    finish do loop finder for  zsidezn0     ====================================================
00784     
00785     //============
00786     
00787     //===========================================================================================================================
00788     
00789     //===========================================================================================================================
00790   }// for zsidezn0 
00791   //===========================================================================================================================
00792   
00793 #ifdef debugsophisticated
00794   std::cout << " numberXtracks= " << numberXtracks << " numberYtracks= " << numberYtracks << std::endl;
00795 #endif
00796   //===========================================================================================================================
00797   //===========================================================================================================================
00798   //===========================================================================================================================
00799 
00800   // case X and Y plane types are available
00801   if(zn0>2) {
00802   //===========================================================================================================================
00803   // match selected X and Y tracks to each other: tgphi=By/Bx->phi=artg(By/Bx); tgtheta=Bx/cosphi=By/sinphi->  ================
00804   //                min of |Bx/cosphi-By/sinphi|                                                               ================
00805 
00806   //  
00807 #ifdef debugsophisticated
00808       std::cout << " numberXtracks= " << numberXtracks << " numberYtracks= " << numberYtracks << std::endl;
00809 #endif
00810       if(numberXtracks>0) {
00811         int newxnum[10], newynum[10];// max # tracks = restracks = 10
00812         int nmathed=0;
00813         do {
00814           double dthmin= 999999.; 
00815           int trminx=-1, trminy=-1;
00816           for (int trx=0; trx<numberXtracks; ++trx) {
00817 #ifdef debugsophisticated
00818             std::cout << "----------- trx= " << trx << " nmathed= " << nmathed << std::endl;
00819 #endif
00820             for (int tr=0; tr<numberYtracks; ++tr) {
00821 #ifdef debugsophisticated
00822               std::cout << "--- tr= " << tr << " nmathed= " << nmathed << std::endl;
00823 #endif
00824               bool YesY=false;
00825               for (int nmx=0; nmx<nmathed; ++nmx) {
00826                 if(trx==newxnum[nmx]) YesY=true;
00827                 if(YesY) break;
00828                 for (int nm=0; nm<nmathed; ++nm) {
00829                   if(tr==newynum[nm]) YesY=true;
00830                   if(YesY) break;
00831                 }
00832               }
00833               if(!YesY) {
00834 //--------------------------------------------------------------------  ----    ----    ----    ----    ----    ----
00835        //double yyyyyy = 999999.;
00836        //if(Bx[trx] != 0.) yyyyyy = Ay[tr]-(Ax[trx]-xxxvtx)*By[tr]/Bx[trx];
00837        //double xxxxxx = 999999.;
00838        //if(By[tr] != 0.) xxxxxx = Ax[trx]-(Ay[tr]-yyyvtx)*Bx[trx]/By[tr];
00839        //double  dthdif= abs(yyyyyy-yyyvtx) + abs(xxxxxx-xxxvtx);
00840 
00841        double  dthdif= abs(AxW[trx]-Ay[tr]) + abs(BxW[trx]-By[tr]);
00842 
00843 #ifdef debugsophisticated
00844        //  std::cout << " yyyyyy= " << yyyyyy << " xxxxxx= " << xxxxxx << " dthdif= " << dthdif << std::endl;
00845   std::cout << " abs(AxW[trx]-Ay[tr]) = " << abs(AxW[trx]-Ay[tr]) << " abs(BxW[trx]-By[tr])= " << abs(BxW[trx]-By[tr]) << " dthdif= " << dthdif << std::endl;
00846 #endif
00847  //--------------------------------------------------------------------     ----        ----    ----    ----    ----    ----
00848                   if( dthdif < dthmin ) {
00849                     dthmin = dthdif;
00850                     trminx = trx;
00851                     trminy = tr;
00852                   }//if  dthdif
00853                   //--------------------------------------------------------------------        
00854               }//if !YesY
00855             }//for y
00856           }// for x
00857           ++nmathed;
00858           if(trminx != -1) {
00859             newxnum[nmathed-1] = trminx;
00860           }
00861           else{
00862             newxnum[nmathed-1] = nmathed-1;
00863           }
00864 #ifdef debugsophisticated
00865           std::cout << " trminx= " << trminx << std::endl;
00866 #endif
00867           if(nmathed>numberYtracks){
00868             newynum[nmathed-1] = -1;
00869 #ifdef debugsophisticated
00870           std::cout << "!!!  nmathed= " << nmathed << " > numberYtracks= " << numberYtracks << std::endl;
00871 #endif
00872           }
00873           else {
00874 #ifdef debugsophisticated
00875             std::cout << " trminy= " << trminy << std::endl;
00876 #endif
00877             newynum[nmathed-1] = trminy;
00878           }    
00879         } while(nmathed<numberXtracks && nmathed < restracks);      
00880         
00881 //
00882 //===========================================================================================================================
00883 //
00884     for (int tr=0; tr<nmathed; ++tr) {
00885       int tx=newxnum[tr];
00886       int ty=newynum[tr];
00887       if(ty==-1){
00888         ty=tx;
00889         Ay[ty]=999.;
00890         By[ty]=999.;
00891         Cy[ty]=999.;
00892         My[ty]=-1;
00893       }//if ty
00894       // test:
00895       //  tx=tr;
00896       //ty=tr;
00897 #ifdef debugsophisticated
00898             if(Mx[tx]>30) {
00899       std::cout << " for track tr= " << tr << " tx= " << tx << " ty= " << ty << std::endl;
00900       std::cout << " Ax= " << Ax[tx]   << " Ay= " << Ay[ty]   << std::endl;
00901       std::cout << " Bx= " << Bx[tx]   << " By= " << By[ty]   << std::endl;
00902       std::cout << " Cx= " << Cx[tx]   << " Cy= " << Cy[ty]   << std::endl;
00903       std::cout << " Mx= " << Mx[tx]   << " My= " << My[ty]   << std::endl;
00904       std::cout << " AxW= " << AxW[tx]   << " AyW= " << AyW[ty]   << std::endl;
00905       std::cout << " BxW= " << BxW[tx]   << " ByW= " << ByW[ty]   << std::endl;
00906       std::cout << " CxW= " << CxW[tx]   << " CyW= " << CyW[ty]   << std::endl;
00907       std::cout << " MxW= " << MxW[tx]   << " MyW= " << MyW[ty]   << std::endl;
00908             }
00909 #endif
00910       //   rhits.push_back( TrackFP420(c0X,c1X,chisqX,nhitplanesY,c0Y,c1Y,chisqY,nhitplanesY) );
00911       rhits.push_back( TrackFP420(Ax[tx],Bx[tx],Cx[tx],Mx[tx],Ay[ty],By[ty],Cy[ty],My[ty]) );
00912     }//for tr
00913     //============================================================================================================
00914   }//in  numberXtracks >0
00915   //============
00916 
00917   }
00918   // case Y plane types are available only
00919   else if(zn0==1) {
00920     for (int ty=0; ty<numberYtracks; ++ty) {
00921 #ifdef debugsophisticated
00922       std::cout << " for track ty= " << ty << std::endl;
00923       std::cout << " Ay= " << Ay[ty]   << std::endl;
00924       std::cout << " By= " << By[ty]   << std::endl;
00925       std::cout << " Cy= " << Cy[ty]   << std::endl;
00926       std::cout << " My= " << My[ty]   << std::endl;
00927       std::cout << " AyW= " << AyW[ty]   << std::endl;
00928       std::cout << " ByW= " << ByW[ty]   << std::endl;
00929       std::cout << " CyW= " << CyW[ty]   << std::endl;
00930       std::cout << " MyW= " << MyW[ty]   << std::endl;
00931 #endif
00932       rhits.push_back( TrackFP420(AyW[ty],ByW[ty],CyW[ty],MyW[ty],Ay[ty],By[ty],Cy[ty],My[ty]) );
00933     }//for ty
00934     //============
00935   }
00936   // case X plane types are available only
00937   else if(zn0==2) {
00938     for (int tx=0; tx<numberXtracks; ++tx) {
00939 #ifdef debugsophisticated
00940       std::cout << " for track tx= " << tx << std::endl;
00941       std::cout << " Ax= " << Ax[tx]   << std::endl;
00942       std::cout << " Bx= " << Bx[tx]   << std::endl;
00943       std::cout << " Cx= " << Cx[tx]   << std::endl;
00944       std::cout << " Mx= " << Mx[tx]   << std::endl;
00945       std::cout << " AxW= " << AxW[tx]   << std::endl;
00946       std::cout << " BxW= " << BxW[tx]   << std::endl;
00947       std::cout << " CxW= " << CxW[tx]   << std::endl;
00948       std::cout << " MxW= " << MxW[tx]   << std::endl;
00949 #endif
00950       rhits.push_back( TrackFP420(Ax[tx],Bx[tx],Cx[tx],Mx[tx],AxW[tx],BxW[tx],CxW[tx],MxW[tx]) );
00951     }//for tx
00952     //============
00953   }//zn0
00954 
00955 
00956 
00957 
00959 
00960 
00961 
00962   return rhits;
00963   //============
00964 }
00968 
00969