Reference Manual

Go to the documentation of this file.
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 arn0, int axytype, double az420, double azD2, double azD3, double apitchX, double apitchY, double apitchXW, double apitchYW, double aZGapLDet, double aZSiStep, double aZSiPlane, double aZSiDet, double azBlade, double agapBlade, bool aUseHalfPitchShiftInX, bool aUseHalfPitchShiftInY, bool aUseHalfPitchShiftInXW, bool aUseHalfPitchShiftInYW, double adXX, double adYY, float achiCutX, float achiCutY, double azinibeg, int verbosity, double aXsensorSize, double aYsensorSize) {
00024   //
00025   // Everything that depend on the det
00026   //
00027   verbos=verbosity;
00028   sn0 = asn0;
00029   pn0 = apn0;
00030   rn0 = arn0;
00031   xytype = axytype;
00032   z420= az420;
00033   zD2 = azD2;
00034   zD3 = azD3;
00035   //zUnit= azUnit;
00036   pitchX = apitchX;
00037   pitchY = apitchY;
00038   pitchXW = apitchXW;
00039   pitchYW = apitchYW;
00040   ZGapLDet = aZGapLDet;
00041   ZSiStep = aZSiStep; 
00042   ZSiPlane = aZSiPlane;
00043   ZSiDet = aZSiDet;
00044   zBlade = azBlade;
00045   gapBlade = agapBlade;
00046 
00047   UseHalfPitchShiftInX = aUseHalfPitchShiftInX;
00048   UseHalfPitchShiftInY = aUseHalfPitchShiftInY;
00049   UseHalfPitchShiftInXW = aUseHalfPitchShiftInXW;
00050   UseHalfPitchShiftInYW = aUseHalfPitchShiftInYW;
00051   dXX = adXX;
00052   dYY = adYY;
00053   chiCutX = achiCutX;
00054   chiCutY = achiCutY;
00055   zinibeg = azinibeg;
00056   XsensorSize = aXsensorSize;
00057   YsensorSize = aYsensorSize;
00058 
00059   if (verbos > 0) {
00060     std::cout << "TrackProducerFP420: call constructor" << std::endl;
00061     std::cout << " sn0= " << sn0 << " pn0= " << pn0 << " rn0= " << rn0 << " xytype= " << xytype << std::endl;
00062     std::cout << " zD2= " << zD2 << " zD3= " << zD3 << " zinibeg= " << zinibeg << std::endl;
00063     //std::cout << " zUnit= " << zUnit << std::endl;
00064     std::cout << " pitchX= " << pitchX << " pitchY= " << pitchY << std::endl;
00065     std::cout << " ZGapLDet= " << ZGapLDet << std::endl;
00066     std::cout << " ZSiStep= " << ZSiStep << " ZSiPlane= " << ZSiPlane << std::endl;
00067     std::cout << " ZSiDet= " <<ZSiDet  << std::endl;
00068     std::cout << " UseHalfPitchShiftInX= " << UseHalfPitchShiftInX << " UseHalfPitchShiftInY= " << UseHalfPitchShiftInY << std::endl;
00069     std::cout << "TrackProducerFP420:----------------------" << std::endl;
00070     std::cout << " dXX= " << dXX << " dYY= " << dYY << std::endl;
00071     std::cout << " chiCutX= " << chiCutX << " chiCutY= " << chiCutY << std::endl;
00072   }
00073 
00074   theFP420NumberingScheme = new FP420NumberingScheme();
00075 
00076 
00077 
00078 }
00080 
00081 
00082 
00083 
00084 
00086 std::vector<TrackFP420> TrackProducerFP420::trackFinderSophisticated(edm::Handle<ClusterCollectionFP420> input, int det){
00088   
00089   std::vector<TrackFP420> rhits;
00090   int restracks = 10;// max # tracks
00091   rhits.reserve(restracks); 
00092   rhits.clear();
00093   double Ax[10]; double Bx[10]; double Cx[10]; int Mx[10];
00094   double Ay[10]; double By[10]; double Cy[10]; int My[10];
00095   double AxW[10]; double BxW[10]; double CxW[10]; int MxW[10];
00096   double AyW[10]; double ByW[10]; double CyW[10]; int MyW[10];
00097   if (verbos > 0) {
00098     std::cout << "===============================================================================" << std::endl; 
00099     std::cout << "=================================================================" << std::endl; 
00100     std::cout << "==========================================================" << std::endl; 
00101     std::cout << "=                                                 =" << std::endl; 
00102     std::cout << "TrackProducerFP420: Start trackFinderSophisticated " << std::endl; 
00103   }
00106 // xytype is the sensor grid arrangment
00107   if( xytype < 1 || xytype > 2 ){
00108     std::cout << "TrackProducerFP420:ERROR in trackFinderSophisticated: check xytype = " << xytype << std::endl; 
00109     return rhits;
00110   }
00111 // sn0= 3 - 2St configuration, sn0= 4 - 3St configuration 
00112 //  if( sn0 < 3 || sn0 > 4 ){
00113   if( sn0 != 3 ){
00114     std::cout << "TrackProducerFP420:ERROR in trackFinderSophisticated: check sn0 (configuration) = " << sn0 << std::endl; 
00115     return rhits;
00116   }
00118   int zbeg = 1, zmax=3;// means layer 1 and 2 in superlayer, i.e. for loop: 1,2
00120   //   .
00121   int reshits1 = 12;// is max # cl in sensor of copyinlayer=1
00122   int reshits2 = 24;// (reshits2-reshits1) is max # cl in sensors of copyinlayers= 2 or 3
00123   //  int resplanes = 20;
00124   int nX[20], nY[20];// resplanes =20 NUMBER OF PLANES; nX, nY - # cl for every X and Y plane
00125   int uX[20], uY[20];// resplanes =20 NUMBER OF PLANES; nX, nY - current # cl used for every X and Y plane
00126   double zX[24][20], xX[24][20], wX[24][20];
00127   double zY[24][20], yY[24][20], wY[24][20];
00128   double             yXW[24][20], wXW[24][20];
00129   double             xYW[24][20], wYW[24][20];
00130   bool qX[24][20], qY[24][20];
00131   //   .
00132   int txf = 0; int txs1 = 0; int txss = 0;
00133   int tyf = 0; int tys1 = 0; int tyss = 0;
00134   //   .
00135   double pitch=0.;
00136   double pitchW=0.;
00137   if(xytype==1){
00138     pitch=pitchY;
00139     pitchW=pitchYW;
00140   }
00141   else if(xytype==2){
00142     pitch=pitchX;
00143     pitchW=pitchXW;
00144   }
00145 
00146 
00147      //current change of geometry:
00148     float Xshift = pitch/2.;
00149     float Yshift = pitchW/2.;
00150     
00151     //
00152     int nmetcurx=0;
00153     int nmetcury=0;
00154     unsigned int ii0 = 999999;
00155     int allplacesforsensors=7;
00156     for (int sector=1; sector < sn0; sector++) {
00157       for (int zmodule=1; zmodule<pn0; zmodule++) {
00158         for (int zsideinorder=1; zsideinorder<allplacesforsensors; zsideinorder++) {
00159           int zside = theFP420NumberingScheme->FP420NumberingScheme::realzside(rn0, zsideinorder);// 1, 3, 5, 2, 4, 6
00160           if (verbos  == -49) {
00161             std::cout << "TrackProducerFP420:  sector= " << sector << " zmodule= " << zmodule << " zsideinorder= " << zsideinorder << " zside= " << zside << " det= " << det << std::endl; 
00162           }
00163           if(zside != 0) {
00164             int justlayer = theFP420NumberingScheme->FP420NumberingScheme::unpackLayerIndex(rn0, zside);// 1, 2
00165             if(justlayer<1||justlayer>2) {
00166               std::cout << "TrackProducerFP420:WRONG  justlayer= " << justlayer << std::endl; 
00167             }
00168             int copyinlayer = theFP420NumberingScheme->FP420NumberingScheme::unpackCopyIndex(rn0, zside);// 1, 2, 3
00169             if(copyinlayer<1||copyinlayer>3) {
00170               std::cout << "TrackProducerFP420:WRONG  copyinlayer= " << copyinlayer << std::endl; 
00171             }
00172             int orientation = theFP420NumberingScheme->FP420NumberingScheme::unpackOrientation(rn0, zside);// Front: = 1; Back: = 2
00173             if(orientation<1||orientation>2) {
00174               std::cout << "TrackProducerFP420:WRONG  orientation= " << orientation << std::endl; 
00175             }
00176             // ii is a continues numbering of planes(!)  over two arm FP420 set up
00177             //                                                                    and  ...[ii] massives have prepared in use of ii
00178             int detfixed=1;// use this treatment for each set up arm, hence no sense to repete the same for +FP420 and -FP420;
00179             int nlayers=3;// 2=3-1
00180             unsigned int ii = theFP420NumberingScheme->FP420NumberingScheme::packMYIndex(nlayers,pn0,sn0,detfixed,justlayer,sector,zmodule)-1;// substruct 1 from 1(+1), 2(+2), 3(+3),4(+4),5...,6...,7...,8...,9...,10... (1st Station)              ,11...,12...,13,...20... (2nd Station)
00181             // ii = 0-19   --> 20 items
00183 
00184             if (verbos == -49) {
00185               std::cout << "TrackProducerFP420:  justlayer= " << justlayer << " copyinlayer= " << copyinlayer << " ii= " << ii << std::endl; 
00186             }
00187             
00188             double zdiststat = 0.;
00189             if(sn0<4) {
00190               if(sector==2) zdiststat = zD3;
00191             }
00192             else {
00193               if(sector==2) zdiststat = zD2;
00194               if(sector==3) zdiststat = zD3;
00195             }
00196             double kplane = -(pn0-1)/2 - 0.5  +  (zmodule-1); //-3.5 +0...5 = -3.5,-2.5,-1.5,+2.5,+1.5
00197             double zcurrent = zinibeg + z420 + (ZSiStep-ZSiPlane)/2  + kplane*ZSiStep + zdiststat;  
00198             //double zcurrent = zinibeg +(ZSiStep-ZSiPlane)/2  + kplane*ZSiStep + (sector-1)*zUnit;  
00199             
00200             if(justlayer==1){
00201              if(orientation==1) zcurrent += (ZGapLDet+ZSiDet/2);
00202              if(orientation==2) zcurrent += zBlade-(ZGapLDet+ZSiDet/2);
00203             }
00204             if(justlayer==2){
00205              if(orientation==1) zcurrent += (ZGapLDet+ZSiDet/2)+zBlade+gapBlade;
00206              if(orientation==2) zcurrent += 2*zBlade+gapBlade-(ZGapLDet+ZSiDet/2);
00207             }
00208             //   .
00209             //
00210             if(det == 2) zcurrent = -zcurrent;
00211             //
00212             //
00213             //   .
00214             // local - global systems with possible shift of every second plate:
00215             
00216             // for xytype=1
00217             float dYYcur = dYY;// XSiDet/2.
00218             float dYYWcur = dXX;//(BoxYshft+dYGap) + (YSi - YSiDet)/2. = 4.7
00219             // for xytype=2
00220             float dXXcur = dXX;//(BoxYshft+dYGap) + (YSi - YSiDet)/2. = 4.7
00221             float dXXWcur = dYY;// XSiDet/2.
00222             //   .
00223             if(justlayer==2) {
00224               // X-type: x-coord
00225               if (UseHalfPitchShiftInX == true){
00226                 dXXcur += Xshift;
00227               }
00228               // X-type: y-coord
00229               if (UseHalfPitchShiftInXW == true){
00230                 dXXWcur -= Yshift;
00231               }
00232             }
00233             //
00234             double XXXDelta = 0.0;
00235             if(copyinlayer==2) { XXXDelta = XsensorSize;}
00236             if(copyinlayer==3) { XXXDelta = 2.*XsensorSize;}
00237             double YYYDelta = 0.0;
00238             if(copyinlayer==2) { YYYDelta = XsensorSize;}
00239             if(copyinlayer==3) { YYYDelta = 2.*XsensorSize;}
00240             //   .
00241             //   GET CLUSTER collection  !!!!
00242             //   .
00243             unsigned int iu=theFP420NumberingScheme->FP420NumberingScheme::packMYIndex(rn0,pn0,sn0,det,zside,sector,zmodule);
00244             if (verbos > 0 ) {
00245               std::cout << "TrackProducerFP420: check        iu = " << iu << std::endl; 
00246               std::cout << "TrackProducerFP420:  sector= " << sector << " zmodule= " << zmodule << " zside= " << zside << " det= " << det << " rn0= " << rn0 << " pn0= " << pn0 << " sn0= " << sn0 << " copyinlayer= " << copyinlayer << std::endl; 
00247             }
00248             //============================================================================================================ put into currentclust
00249             std::vector<ClusterFP420> currentclust;
00250             currentclust.clear();
00251             ClusterCollectionFP420::Range outputRange;
00252             outputRange = input->get(iu);
00253             // fill output in currentclust vector (for may be sorting? or other checks)
00254             ClusterCollectionFP420::ContainerIterator sort_begin = outputRange.first;
00255             ClusterCollectionFP420::ContainerIterator sort_end = outputRange.second;
00256             for ( ;sort_begin != sort_end; ++sort_begin ) {
00257               //  std::sort(currentclust.begin(),currentclust.end());
00258               currentclust.push_back(*sort_begin);
00259             } // for
00260             
00261             if (verbos > 0 ) {
00262               std::cout << "TrackProducerFP420: currentclust.size = " << currentclust.size() << std::endl; 
00263             }
00264             //============================================================================================================
00265             
00266             std::vector<ClusterFP420>::const_iterator simHitIter = currentclust.begin();
00267             std::vector<ClusterFP420>::const_iterator simHitIterEnd = currentclust.end();
00268             
00269             if(xytype ==1){
00270               if(ii != ii0) {
00271                 ii0=ii;
00272                 nY[ii] = 0;// # cl in every Y plane (max is reshits)
00273                 uY[ii] = 0;// current used # cl in every X plane 
00274                 nmetcury=0;
00275               }
00276             }
00277             else if(xytype ==2){
00278               if(ii != ii0) {
00279                 ii0=ii;
00280                 nX[ii] = 0;// # cl in every X plane (max is reshits)
00281                 uX[ii] = 0;// current used # cl in every X plane 
00282                 nmetcurx=0;
00283               }
00284             }
00285             // loop in #clusters of current sensor
00286             for (;simHitIter != simHitIterEnd; ++simHitIter) {
00287               const ClusterFP420 icluster = *simHitIter;
00288               
00289               // fill vectors for track reconstruction
00290               
00291               //disentangle complicated pattern recognition of hits?
00292               // Y:
00293               if(xytype ==1){
00294                 nY[ii]++;               
00295                 if(copyinlayer==1 && nY[ii]>reshits1){
00296                   nY[ii]=reshits1;
00297                   std::cout << "WARNING-ERROR:TrackproducerFP420: currentclust.size()= " << currentclust.size() <<" bigger reservated number of hits" << " zcurrent=" << zY[nY[ii]-1][ii] << " copyinlayer= "  << copyinlayer << " ii= "  << ii << std::endl;
00298                 }
00299                 if(copyinlayer !=1 && nY[ii]>reshits2){
00300                   nY[ii]=reshits2;
00301                   std::cout << "WARNING-ERROR:TrackproducerFP420: currentclust.size()= " << currentclust.size() <<" bigger reservated number of hits" << " zcurrent=" << zY[nY[ii]-1][ii] << " copyinlayer= "  << copyinlayer << " ii= "  << ii << std::endl;
00302                 }
00303                 zY[nY[ii]-1][ii] = zcurrent;
00304                 yY[nY[ii]-1][ii] = icluster.barycenter()*pitch+0.5*pitch+YYYDelta;
00305                 xYW[nY[ii]-1][ii] = icluster.barycenterW()*pitchW+0.5*pitchW;
00306                 // go to global system:
00307                 yY[nY[ii]-1][ii] = yY[nY[ii]-1][ii] - dYYcur; 
00308                 wY[nY[ii]-1][ii] = 1./(icluster.barycerror()*pitch);//reciprocal of the variance for each datapoint in y
00309                 wY[nY[ii]-1][ii] *= wY[nY[ii]-1][ii];//reciprocal of the variance for each datapoint in y
00310                 if(det == 2) {
00311                   xYW[nY[ii]-1][ii] =(xYW[nY[ii]-1][ii]+dYYWcur); 
00312                 }
00313                 else {
00314                   xYW[nY[ii]-1][ii] =-(xYW[nY[ii]-1][ii]+dYYWcur); 
00315                 }
00316                 wYW[nY[ii]-1][ii] = 1./(icluster.barycerrorW()*pitchW);//reciprocal of the variance for each datapoint in y
00317                 wYW[nY[ii]-1][ii] *= wYW[nY[ii]-1][ii];//reciprocal of the variance for each datapoint in y
00318                 qY[nY[ii]-1][ii] = true;
00319                 if(copyinlayer==1 && nY[ii]==reshits1) break;
00320                 if(copyinlayer !=1 && nY[ii]==reshits2) break;
00321               }
00322               // X:
00323               else if(xytype ==2){
00324                 nX[ii]++;       
00325                 if (verbos == -49) {
00326                   std::cout << "TrackproducerFP420: nX[ii]= " << nX[ii] << " Ncl= " << currentclust.size() << " copyinlayer= "  << copyinlayer << " ii= " << ii << " zcurrent = " << zcurrent << " xX= " << icluster.barycenter()*pitch+0.5*pitch+XXXDelta << " yXW= " << icluster.barycenterW()*pitchW+0.5*pitchW << " det= " << det << " cl.size= " << icluster.amplitudes().size() << " cl.ampl[0]= " << icluster.amplitudes()[0] << std::endl;
00327                 }
00328                 if(copyinlayer==1 && nX[ii]>reshits1){
00329                   std::cout << "WARNING-ERROR:TrackproducerFP420: nX[ii]= " << nX[ii] <<" bigger reservated number of hits" << " currentclust.size()= " << currentclust.size() << " copyinlayer= "  << copyinlayer << " ii= " << ii << std::endl;
00330                   nX[ii]=reshits1;
00331                 }
00332                 if(copyinlayer !=1 && nX[ii]>reshits2){
00333                   std::cout << "WARNING-ERROR:TrackproducerFP420: nX[ii]= " << nX[ii] <<" bigger reservated number of hits" << " currentclust.size()= " << currentclust.size() << " copyinlayer= "  << copyinlayer << " ii= " << ii << std::endl;
00334                   nX[ii]=reshits2;
00335                 }
00336                 zX[nX[ii]-1][ii] = zcurrent;
00337                 xX[nX[ii]-1][ii] = icluster.barycenter()*pitch+0.5*pitch+XXXDelta;
00338                 yXW[nX[ii]-1][ii] = icluster.barycenterW()*pitchW+0.5*pitchW;
00339                 // go to global system:
00340                 xX[nX[ii]-1][ii] =-(xX[nX[ii]-1][ii]+dXXcur); 
00341                 wX[nX[ii]-1][ii] = 1./(icluster.barycerror()*pitch);//reciprocal of the variance for each datapoint in y
00342                 wX[nX[ii]-1][ii] *= wX[nX[ii]-1][ii];//reciprocal of the variance for each datapoint in y
00343                 if(det == 2) {
00344                   yXW[nX[ii]-1][ii] = -(yXW[nX[ii]-1][ii] - dXXWcur); 
00345                 }
00346                 else {
00347                   yXW[nX[ii]-1][ii] = yXW[nX[ii]-1][ii] - dXXWcur; 
00348                 }
00349                 wXW[nX[ii]-1][ii] = 1./(icluster.barycerrorW()*pitchW);//reciprocal of the variance for each datapoint in y
00350                 wXW[nX[ii]-1][ii] *= wXW[nX[ii]-1][ii];//reciprocal of the variance for each datapoint in y
00351                 qX[nX[ii]-1][ii] = true;
00352                 if (verbos == -29) {
00353                   std::cout << "trackFinderSophisticated: nX[ii]= " << nX[ii]<< " ii = " << ii << " zcurrent = " << zcurrent << " yXW[nX[ii]-1][ii] = " << yXW[nX[ii]-1][ii] << " xX[nX[ii]-1][ii] = " << xX[nX[ii]-1][ii] << std::endl;
00354                   std::cout << "  XXXDelta= " << XXXDelta << "  dXXcur= " << dXXcur << "  -dXXWcur= " << -dXXWcur << std::endl;
00355                   std::cout << "  icluster.barycerrorW()*pitchW= " << icluster.barycerrorW()*pitchW << "  wXW[nX[ii]-1][ii]= " <<wXW[nX[ii]-1][ii]  << std::endl;
00356                   std::cout << " -icluster.barycenterW()*pitchW+0.5*pitchW = " << icluster.barycenterW()*pitchW+0.5*pitchW << std::endl;
00357                   std::cout << "============================================================" << std::endl;
00358                 }
00359                 if (verbos  > 0) {
00360                   std::cout << "trackFinderSophisticated: nX[ii]= " << nX[ii]<< " ii = " << ii << " zcurrent = " << zcurrent << " xX[nX[ii]-1][ii] = " << xX[nX[ii]-1][ii] << std::endl;
00361                   std::cout << " wX[nX[ii]-1][ii] = " << wX[nX[ii]-1][ii] << " wXW[nX[ii]-1][ii] = " << wXW[nX[ii]-1][ii] << std::endl;
00362                   std::cout << " -icluster.barycenter()*pitch-0.5*pitch = " << -icluster.barycenter()*pitch-0.5*pitch << " -dXXcur = " << -dXXcur << " -XXXDelta = " << -XXXDelta << std::endl;
00363                   std::cout << "============================================================" << std::endl;
00364                 }
00365 
00366                 if(copyinlayer==1 && nX[ii]==reshits1) break;
00367                 if(copyinlayer !=1 && nX[ii]==reshits2) break;
00368               }// if(xytype
00369               
00370             } // for loop in #clusters (can be breaked)
00371             
00372             // Y:
00373             if(xytype ==1){
00374               if(nY[ii] > nmetcury) {  /* # Y-planes w/ clusters */
00375                 nmetcury=nY[ii];
00376                 ++tyf; if(sector==1) ++tys1; if(sector==(sn0-1)) ++tyss;
00377               }   
00378             }
00379             // X:
00380             else if(xytype ==2){
00381               if(nX[ii] > nmetcurx) {  /* # X-planes w/ clusters */
00382                 nmetcurx=nX[ii];
00383                 ++txf; if(sector==1) ++txs1; if(sector==(sn0-1)) ++txss;
00384               }   
00385             }
00386             //================================== end of for loops in continuius number iu:
00387           }//if(zside!=0
00388         }   // for zsideinorder
00389       }   // for zmodule
00390     }   // for sector
00391     if (verbos > 0) {
00392       std::cout << "trackFinderSophisticated: tyf= " << tyf<< " tys1 = " << tys1 << " tyss = " << tyss << std::endl;
00393       std::cout << "trackFinderSophisticated: txf= " << txf<< " txs1 = " << txs1 << " txss = " << txss << std::endl;
00394       std::cout << "============================================================" << std::endl;
00395     }
00396     
00397     //===========================================================================================================================
00398     //===========================================================================================================================
00399     //===========================================================================================================================
00400     //======================    start road finder   =============================================================================
00401     //===========================================================================================================================
00402 
00403   //  int nitMax=5;// max # iterations to find track
00404   int nitMax=10;// max # iterations to find track using go over of different XZ and YZ fits to find the good chi2X and chi2Y simultaneously(!!!) 
00405 
00406   // criteria for track selection: 
00407   // track is selected if for 1st station #cl >=pys1Cut
00408   //  int  pys1Cut = 5, pyssCut = 5, pyallCut=12;
00409   //  int  pys1Cut = 1, pyssCut = 1, pyallCut= 3;
00410 
00411 //  int  pys1Cut = 3, pyssCut = 3, pyallCut= 6; // before geom. update
00412 //  int  pys1Cut = 2, pyssCut = 2, pyallCut= 4; // bad for 5 layers per station
00413   int  pys1Cut = 3, pyssCut = 1, pyallCut= 5;
00414   
00415   //  double yyyvtx = 0.0, xxxvtx = -15;  //mm
00416   
00418   //
00419   // for configuration: 3St, 1m for 1-2 St:
00420   // double sigman=0.1, ssigma = 1.0, sigmam=0.15;/* ssigma is foreseen to match 1st point of 2nd Station*/
00421   //
00422   // for equidistant 3 Stations:
00423   //
00424   // for tests:
00425   //  double sigman=118., ssigma = 299., sigmam=118.;
00426   // RMS1=0.013, RMS2 = 1.0, RMS3 = 0.018 see plots d1XCL, d2XCL, d3XCL
00427   //
00428   //  double sigman=0.05, ssigma = 2.5, sigmam=0.06;
00429   //  double sigman=0.18, ssigma = 1.8, sigmam=0.18;
00430   //  double sigman=0.18, ssigma = 2.9, sigmam=0.18;
00431   //
00433   // for 3 Stations:
00434   // LAST:
00435   double sigman=0.18, ssigma = 2.5, sigmam=0.18;
00436   if( sn0 < 4 ){
00437     // for 2 Stations:
00438     // sigman=0.24, ssigma = 4.2, sigmam=0.33;
00439     //  sigman=0.18, ssigma = 3.9, sigmam=0.18;
00440     // sigman=0.18, ssigma = 3.6, sigmam=0.18;
00441     //
00442 
00443     //
00444 
00445     //  sigman=0.18, ssigma = 3.3, sigmam=0.18;// before geometry update for 4 sensors per superlayer
00446     //    sigman=0.30, ssigma = 7.1, sigmam=0.40;// for update
00447      sigman=0.30, ssigma = 8.0, sigmam=1.0;// for matching update to find point nearby to fit track in 1st plane of 2nd Station 
00448     //
00449   }
00450   if (verbos > 0) {
00451     std::cout << "trackFinderSophisticated: ssigma= " << ssigma << std::endl;
00452   }
00453   
00455   
00456   /* ssigma = 3. * 8000.*(0.025+0.009)/sqrt(pn0-1)/100. = 2.9 mm(!!!)----
00457      ssigma is reduced by factor k_reduced = (sn0-1)-sector+1 = sn0-sector
00458      # Stations  currentStation
00459      2Stations:     sector=2,         sn0=3 , sn0-sector = 1 --> k_reduced = 1
00460      3Stations:     sector=2,         sn0=4 , sn0-sector = 2 --> k_reduced = 2
00461      3Stations:     sector=3,         sn0=4 , sn0-sector = 1 --> k_reduced = 1
00462   */
00463   int numberXtracks=0, numberYtracks=0, totpl = 2*(pn0-1)*(sn0-1); double sigma;
00464 
00465   for (int xytypecurrent=xytype; xytypecurrent<xytype+1; ++xytypecurrent) {
00466     if (verbos > 0) {
00467       std::cout << "trackFinderSophisticated: xytypecurrent= " << xytypecurrent << std::endl;
00468     }
00469     
00470     //
00471     //
00472     double tg0 = 0.;
00473     int qAcl[20], qAii[20], fip=0, niteration = 0;
00474     int ry = 0, rys1 = 0, ryss = 0;
00475     int tas1=tys1, tass=tyss, taf=tyf;
00476     bool SelectTracks = true;
00477     //
00479   //   .
00480 
00481   double yA[24][20], zA[24][20], wA[24][20]; int nA[20], uA[20]; bool qA[24][20];
00482     //
00483     // Y:
00484   //======================    start road finder  for xytypecurrent = 1      ===========================================================
00485     if(xytypecurrent ==1){
00486   //===========================================================================================================================
00487       numberYtracks=0;  
00488       tg0= 3*1./(800.+20.); // for Y: 1cm/...   *3 - 3sigma range
00489       tas1=tys1;
00490       tass=tyss;
00491       taf=tyf;
00492       for (int ii=0; ii < totpl; ++ii) {
00493         if (verbos > 0) {
00494           std::cout << "trackFinderSophisticated: ii= " << ii << " nY[ii]= " << nY[ii] << std::endl;
00495           std::cout << "trackFinderSophisticated: ii= " << ii << " uY[ii]= " << uY[ii] << std::endl;
00496         }
00497         nA[ii] = nY[ii];
00498         uA[ii] = uY[ii];
00499         for (int cl=0; cl<nA[ii]; ++cl) {
00500           if (verbos > 0) {
00501             std::cout << " cl= " << cl << " yY[cl][ii]= " << yY[cl][ii] << std::endl;
00502             std::cout << " zY[cl][ii]= " << zY[cl][ii] << " wY[cl][ii]= " << wY[cl][ii] << " qY[cl][ii]= " << qY[cl][ii] << std::endl;
00503           }
00504           yA[cl][ii] = yY[cl][ii];
00505           zA[cl][ii] = zY[cl][ii];
00506           wA[cl][ii] = wY[cl][ii];
00507           qA[cl][ii] = qY[cl][ii];
00508         }
00509       }
00510   //===========================================================================================================================
00511     }// if xytypecurrent ==1
00512     // X:
00513   //======================    start road finder  for superlayer = 2      ===========================================================
00514     else if(xytypecurrent ==2){
00515   //===========================================================================================================================
00516       numberXtracks=0;  
00517       tg0= 3*2./(800.+20.); // for X: 2cm/...   *3 - 3sigma range
00518       tas1=txs1;
00519       tass=txss;
00520       taf=txf;
00521       for (int ii=0; ii < totpl; ++ii) {
00522         if (verbos > 0) {
00523           std::cout << "trackFinderSophisticated: ii= " << ii << " nX[ii]= " << nX[ii] << std::endl;
00524           std::cout << "trackFinderSophisticated: ii= " << ii << " uX[ii]= " << uX[ii] << std::endl;
00525         }
00526         nA[ii] = nX[ii];
00527         uA[ii] = uX[ii];
00528         for (int cl=0; cl<nA[ii]; ++cl) {
00529           if (verbos == -29) {
00530             std::cout << " cl= " << cl << " xX[cl][ii]= " << xX[cl][ii] << std::endl;
00531             std::cout << " zX[cl][ii]= " << zX[cl][ii] << " wX[cl][ii]= " << wX[cl][ii] << " qX[cl][ii]= " << qX[cl][ii] << std::endl;
00532           }
00533           yA[cl][ii] = xX[cl][ii];
00534           zA[cl][ii] = zX[cl][ii];
00535           wA[cl][ii] = wX[cl][ii];
00536           qA[cl][ii] = qX[cl][ii];
00537         }
00538       }
00539   //===========================================================================================================================
00540     }// if xytypecurrent ==xytype
00541 
00542 
00543     
00544   //======================    start road finder        ====================================================
00545     if (verbos > 0) {
00546       std::cout << "                  start road finder                        " << std::endl;
00547     }
00548     do {
00549       double fyY[20], fzY[20], fwY[20];
00550       double fyYW[20],         fwYW[20];
00551       int py = 0, pys1 = 0, pyss = 0;
00552       bool NewStation = false, py1first = false;
00553       for (int sector=1; sector < sn0; ++sector) {
00554         double tav=0., t1=0., t2=0., t=0., sm;
00555         int stattimes=0;
00556         if( sector != 1 ) {
00557           NewStation = true;  
00558         }
00559         for (int zmodule=1; zmodule<pn0; ++zmodule) {
00560           for (int justlayer=zbeg; justlayer<zmax; justlayer++) {
00561             // iu is a continues numbering of planes(!) 
00562             int detfixed=1;// use this treatment for each set up arm, hence no sense to do it differently for +FP420 and -FP420;
00563             int nlayers=3;// 2=3-1
00564             unsigned int ii = theFP420NumberingScheme->FP420NumberingScheme::packMYIndex(nlayers,pn0,sn0,detfixed,justlayer,sector,zmodule)-1;// substruct 1 from 1(+1), 2(+2), 3(+3),4(+4),5...,6...,7...,8...,9...,10... (1st Station)              ,11...,12...,13,...20... (2nd Station)
00565             // ii = 0-19   --> 20 items
00566             
00567             if(nA[ii]!=0  && uA[ii]!= nA[ii]) { 
00568               
00569               ++py; if(sector==1) ++pys1; if(sector==(sn0-1)) ++pyss;
00570               if(py==2 && sector==1) { 
00571                 // find closest cluster in X                   .
00572                 double dymin=9999999., df2; int cl2=-1;
00573                 for (int cl=0; cl<nA[ii]; ++cl) {
00574                   if(qA[cl][ii]){
00575                     df2 = std::abs(fyY[fip]-yA[cl][ii]);
00576                     if(df2 < dymin) {
00577                       dymin = df2;
00578                       cl2=cl;
00579                     }//if(df2           
00580                   }//if(qA              
00581                 }//for(cl
00582                 // END of finding of closest cluster in X                   .
00583                 if(cl2!=-1){
00584                   t=(yA[cl2][ii]-fyY[fip])/(zA[cl2][ii]-fzY[fip]);
00585                   t1 = t*wA[cl2][ii];
00586                   t2 = wA[cl2][ii];
00587                   if (verbos > 0) {
00588                     std::cout << " t= " << t << " tg0= " << tg0 << std::endl;
00589                   }
00590                   if(std::abs(t)<tg0) { 
00591                     qA[cl2][ii] = false;//point is taken, mark it for not using again
00592                     fyY[py-1]=yA[cl2][ii];
00593                     fzY[py-1]=zA[cl2][ii];
00594                     fwY[py-1]=wA[cl2][ii];
00595                     qAcl[py-1] = cl2;
00596                     qAii[py-1] = ii;
00597                     ++uA[ii];
00598                     if (verbos > 0) {
00599                       std::cout << " point is taken, mark it for not using again uA[ii]= " << uA[ii] << std::endl;
00600                     }
00601                     if(uA[ii]==nA[ii]){/* no points anymore for this plane */
00602                       ++ry; if(sector==1) ++rys1; if(sector==(sn0-1)) ++ryss;
00603                     }//if(uA
00604                   }//if abs
00605                   else{
00606                     py--; if(sector==1) pys1--;  if(sector==(sn0-1)) pyss--;
00607                     t1 -= t*wA[cl2][ii]; t2 -= wA[cl2][ii];
00608                   }//if(abs
00609                 }//if(cl2!=-1
00610                 else{
00611                   py--; if(sector==1) pys1--;  if(sector==(sn0-1)) pyss--;
00612                 }//if(cl2!=-1
00613               }//if(py==2
00614               else {
00615                 //                                                                                                                 .
00616                 bool clLoopTrue = true;
00617                 int clcurr=-1;
00618                 for (int clind=0; clind<nA[ii]; ++clind) {
00619                   if(clLoopTrue) {
00620                     int cl=clind;
00621                     if(qA[cl][ii]){
00622                       clcurr = cl;
00623                       if(py<3 ){
00624                         if(py==1) { 
00625                           py1first = true;
00626                           fip=py-1;
00627                           qA[cl][ii] = false;//point is taken, mark it for not using again
00628                           fyY[py-1]=yA[cl][ii];
00629                           fzY[py-1]=zA[cl][ii];
00630                           fwY[py-1]=wA[cl][ii];
00631                           qAcl[py-1] = cl;
00632                           qAii[py-1] = ii;
00633                           ++uA[ii];
00634                           if (verbos > 0) std::cout << " point is taken, mark it uA[ii]= " << uA[ii] << std::endl;
00635                         }//if py=1
00636                         if(uA[ii]==nA[ii]){/* no points anymore for this plane */
00637                           ++ry; if(sector==1) ++rys1; if(sector==(sn0-1)) ++ryss;
00638                         }//if(uA
00639                       }//py<3
00640                       else {
00641                         if(NewStation){
00642                           if( sn0 < 4 ) {
00643                             // stattimes=0 case (point of 1st plane to be matched in new Station)
00644                             sigma = ssigma;
00645                           }
00646                           else {
00647                             sigma = ssigma/(sn0-1-sector);
00648                           }
00649                           //sigma = ssigma/(sn0-sector);
00650                           //if(stattimes==1 || sector==3 ) sigma = msigma * sqrt(1./wA[cl][ii]);
00651                           if(stattimes==1 || sector==3 ) sigma = sigmam; // (1st $3rd Stations for 3St. configur. ), 1st only for 2St. conf.
00652                           //    if(stattimes==1 || sector==(sn0-1) ) sigma = sigmam;
00653                           
00654                           double cov00, cov01, cov11, c0Y, c1Y, chisqY;
00655                           gsl_fit_wlinear (fzY, 1, fwY, 1, fyY, 1, py-1, 
00656                                            &c0Y, &c1Y, &cov00, &cov01, &cov11, 
00657                                            &chisqY);
00658                           
00659                           // find closest cluster in X                   .
00660                           int cl2match=-1;
00661                           double dymin=9999999., df2; 
00662                           for (int clmatch=0; clmatch<nA[ii]; ++clmatch) {
00663                             if(qA[clmatch][ii]){
00664                               double smmatch = c0Y+ c1Y*zA[clmatch][ii];
00665                               df2 = std::abs(smmatch-yA[clmatch][ii]);
00666                               if(df2 < dymin) {
00667                                 dymin = df2;
00668                                 cl2match=clmatch;
00669                               }//if(df2         
00670                             }//if(qA            
00671                           }//for(clmatch
00672                           
00673                           if(cl2match != -1) {
00674                             cl=cl2match;
00675                             clLoopTrue = false; // just not continue the clinid loop
00676                           }
00677                           
00678                           sm = c0Y+ c1Y*zA[cl][ii];
00679                           
00680                           if (verbos > 0) {
00681                             std::cout << " sector= " << sector << " sn0= " << sn0 << " sigma= " << sigma << std::endl;
00682                             std::cout << " stattimes= " << stattimes << " ssigma= " << ssigma << " sigmam= " << sigmam << std::endl;
00683                             std::cout << " sm= " << sm << " c0Y= " << c0Y << " c1Y= " << c1Y << " chisqY= " << chisqY << std::endl;
00684                             std::cout << " zA[cl][ii]= " << zA[cl][ii] << " ii= " << ii << " cl= " << cl << std::endl;
00685                             for (int ct=0; ct<py-1; ++ct) {
00686                               std::cout << " py-1= " << py-1 << " fzY[ct]= " << fzY[ct] << std::endl;
00687                               std::cout << " fyY[ct]= " << fyY[ct] << " fwY[ct]= " << fwY[ct] << std::endl;
00688                             }
00689                           }
00690                           
00691                         }//NewStation 1
00692                         else{
00693                           t=(yA[cl][ii]-fyY[fip])/(zA[cl][ii]-fzY[fip]);
00694                           t1 += t*wA[cl][ii];
00695                           t2 += wA[cl][ii];
00696                           tav=t1/t2;
00697                           sm = fyY[fip]+tav*(zA[cl][ii]-fzY[fip]);
00698                           //sigma = nsigma * sqrt(1./wA[cl][ii]);
00699                           sigma = sigman;
00700                         }
00701                         
00702                         double diffpo = yA[cl][ii]-sm;
00703                         if (verbos > 0) {
00704                           std::cout << " diffpo= " << diffpo << " yA[cl][ii]= " << yA[cl][ii] << " sm= " << sm << " sigma= " << sigma << std::endl;
00705                         }
00706                         
00707                         if(std::abs(diffpo) < sigma ) {
00708                           if(NewStation){
00709                             ++stattimes;
00710                             if(stattimes==1) {
00711                               fip=py-1;
00712                               t1 = 0; t2 = 0;
00713                             }
00714                             else if(stattimes==2){
00715                               NewStation = false; 
00716                               t=(yA[cl][ii]-fyY[fip])/(zA[cl][ii]-fzY[fip]);
00717                               //t1 += t*wA[cl][ii];
00718                               //t2 += wA[cl][ii];
00719                               t1 = t*wA[cl][ii];
00720                               t2 = wA[cl][ii];
00721                             }//if(stattime
00722                           }//if(NewStation 2
00723                           fyY[py-1]=yA[cl][ii];
00724                           fzY[py-1]=zA[cl][ii];
00725                           fwY[py-1]=wA[cl][ii];
00726                           qA[cl][ii] = false;//point is taken, mark it for not using again
00727                           qAcl[py-1] = cl;
00728                           qAii[py-1] = ii;
00729                           ++uA[ii];
00730                           if (verbos > 0) {
00731                             std::cout << " 3333 point is taken, mark it uA[ii]= " << uA[ii] << std::endl;
00732                           }
00733                           if(uA[ii]==nA[ii]){/* no points anymore for this plane */
00734                             ++ry; if(sector==1) ++rys1; if(sector==(sn0-1)) ++ryss;
00735                           }//if(cl==
00736                           //  break; // to go on next plane
00737                         }//if abs
00738                         else{
00739                           t1 -= t*wA[cl][ii]; t2 -= wA[cl][ii];
00740                         }//if abs
00741                       }// if py<3 and else py>3
00742                       
00743                       if(!qA[cl][ii]) break;// go on next plane if point is found among clusters of current plane;
00744                     }// if qA
00745                   } // if clLoopTrue
00746                 }// for cl     --  can be break and return to "for zmodule"
00747                 //                                                                                                                 .
00748                 if( (py!=1 && clcurr != -1 && qA[clcurr][ii]) || (py==1 && !py1first)) { 
00749                   // if point is not found - continue natural loop, but reduce py 
00750                   py--; if(sector==1) pys1--;  if(sector==(sn0-1)) pyss--;
00751                 }//if(py!=1
00752               }//if(py==2 else 
00753             }//if(nA !=0           : inside  this if( -  ask  ++py
00754           }// for justlayer
00755         }// for zmodule
00756       }// for sector
00757       //============
00759       
00760       if (verbos > 0) {
00761         std::cout << "END: pys1= " << pys1 << " pyss = " << pyss << " py = " << py << std::endl;
00762       }
00763       // apply criteria for track selection: 
00764       // do not take track if 
00765       if( pys1 < pys1Cut || pyss < pyssCut || py < pyallCut ){
00766         //      if( pys1<3 || pyss<2 || py<4 ){
00767       }
00768       // do fit:
00769       else{
00771         double cov00, cov01, cov11;
00772         double c0Y, c1Y, chisqY;
00773         gsl_fit_wlinear (fzY, 1, fwY, 1, fyY, 1, py, 
00774                          &c0Y, &c1Y, &cov00, &cov01, &cov11, 
00775                          &chisqY);
00777         // collect cases where the nearby points with the same coordinate exists
00778 //      int pyrepete=py+2;
00779 //      if(py < 11 && chisqY/(py-2) < 0.5) {
00780 //        double fyYold=999999.;
00781 //        for (int ipy=0; ipy<py; ++ipy) {
00782 //          if( fyY[ipy]!=fyYold) --pyrepete;
00783 //          fyYold=fyY[ipy];
00784 //        }
00785 //      }
00787         float chindfx;
00788         if(py>2) {
00789           chindfx = chisqY/(py-2);
00790         }
00791         else{
00792           //      chindfy = chisqY;
00793           chindfx = 9999;
00794         }//py
00795         if (verbos  > 0) {
00796           //      std::cout << " Do FIT XZ: chindfx= " << chindfx << " chisqY= " << chisqY << " py= " << py << " pyrepete= " << pyrepete << std::endl;
00797           std::cout << " Do FIT XZ: chindfx= " << chindfx << " chisqY= " << chisqY << " py= " << py << std::endl;
00798         }
00799         
00801         if (verbos > 0) {
00802           std::cout << " preparation for second order fit for Wide pixels= " << std::endl;
00803         }
00804         for (int ipy=0; ipy<py; ++ipy) {
00805           if(xytypecurrent ==1){
00806             fyYW[ipy]=xYW[qAcl[ipy]][qAii[ipy]];
00807             fwYW[ipy]=wYW[qAcl[ipy]][qAii[ipy]];
00808             if (verbos > 0) {
00809               std::cout << " ipy= " << ipy << std::endl;
00810               std::cout << " qAcl[ipy]= " << qAcl[ipy] << " qAii[ipy]= " << qAii[ipy] << std::endl;
00811               std::cout << " fyYW[ipy]= " << fyYW[ipy] << " fwYW[ipy]= " << fwYW[ipy] << std::endl;
00812             }
00813           }
00814           else if(xytypecurrent ==2){
00815             fyYW[ipy]=yXW[qAcl[ipy]][qAii[ipy]];
00816             fwYW[ipy]=wXW[qAcl[ipy]][qAii[ipy]];
00817             if (verbos ==-29) {
00818               std::cout << " ipy= " << ipy << std::endl;
00819               std::cout << " qAcl[ipy]= " << qAcl[ipy] << " qAii[ipy]= " << qAii[ipy] << std::endl;
00820               std::cout << " fyYW[ipy]= " << fyYW[ipy] << " fwYW[ipy]= " << fwYW[ipy] << std::endl;
00821             }
00822           }
00823         }// for
00824 
00825 
00826 
00827         if (verbos > 0) {
00828           std::cout << " start second order fit for Wide pixels= " << std::endl;
00829         }
00830         double wov00, wov01, wov11;
00831         double w0Y, w1Y, whisqY;
00832         gsl_fit_wlinear (fzY, 1, fwYW, 1, fyYW, 1, py, 
00833                          &w0Y, &w1Y, &wov00, &wov01, &wov11, 
00834                          &whisqY);
00836 
00837 
00838 
00839         float chindfy;
00840         if(py>2) {
00841           chindfy = whisqY/(py-2);
00842         }
00843         else{
00844           //      chindfy = chisqY;
00845           chindfy = 9999;
00846         }//py
00847         
00848         if (verbos > 0) {
00849           std::cout << " chindfy= " << chindfy << " chiCutY= " << chiCutY << std::endl;
00850         }
00851 
00852 
00853         if(xytypecurrent ==1){
00854           if(chindfx < chiCutX && chindfy < chiCutY) {
00855             ++numberYtracks;
00856             Ay[numberYtracks-1] = c0Y; 
00857             By[numberYtracks-1] = c1Y; 
00858             Cy[numberYtracks-1] = chisqY; 
00859             //  My[numberYtracks-1] = py-pyrepete;
00860             My[numberYtracks-1] = py;
00861             AyW[numberYtracks-1] = w0Y; 
00862             ByW[numberYtracks-1] = w1Y; 
00863             CyW[numberYtracks-1] = whisqY; 
00864             MyW[numberYtracks-1] = py;
00865             if (verbos > 0) {
00866               if(py>20) {
00867                 std::cout << " niteration = " << niteration << std::endl;
00868                 std::cout << " chindfy= " << chindfy << " py= " << py << std::endl;
00869                 std::cout << " c0Y= " << c0Y << " c1Y= " << c1Y << std::endl;
00870                 std::cout << " pys1= " << pys1 << " pyss = " << pyss << std::endl;
00871               }
00872             }
00873           }//chindfy
00874         }
00875         else if(xytypecurrent ==2){
00876           if(chindfx < chiCutX && chindfy < chiCutY) {
00877             ++numberXtracks;
00878             Ax[numberXtracks-1] = c0Y; 
00879             Bx[numberXtracks-1] = c1Y; 
00880             Cx[numberXtracks-1] = chisqY; 
00881             //   Mx[numberXtracks-1] = py-pyrepete;
00882             Mx[numberXtracks-1] = py;
00883             AxW[numberXtracks-1] = w0Y; 
00884             BxW[numberXtracks-1] = w1Y; 
00885             CxW[numberXtracks-1] = whisqY; 
00886             MxW[numberXtracks-1] = py;
00887             if (verbos > 0) {
00888               std::cout << " niteration = " << niteration << std::endl;
00889               std::cout << " chindfx= " << chindfy << " px= " << py << std::endl;
00890               std::cout << " c0X= " << c0Y << " c1X= " << c1Y << std::endl;
00891               std::cout << " pxs1= " << pys1 << " pxss = " << pyss << std::endl;
00892             }
00893           }//chindfy
00894         }
00895         
00896         
00897       }//  if else
00899       // do not select tracks anymore if
00900       if (verbos  > 0) {
00901         std::cout << "Current iteration, niteration >= " << niteration << std::endl;
00902         std::cout << " numberYtracks= " << numberYtracks << std::endl;
00903         std::cout << " numberXtracks= " << numberXtracks << std::endl;
00904         std::cout << " pys1= " << pys1 << " pyss = " << pyss << " py = " << py << std::endl;
00905         std::cout << " tas1= " << tas1 << " tass = " << tass << " taf = " << taf << std::endl;
00906         std::cout << " rys1= " << rys1 << " ryss = " << ryss << " ry = " << ry << std::endl;
00907         std::cout << " tas1-rys1= " << tas1-rys1 << " tass-ryss = " << tass-ryss << " taf-ry = " << taf-ry << std::endl;
00908         std::cout << "---------------------------------------------------------- " << std::endl;
00909       }
00910       // let's decide: do we continue track finder procedure
00911       if( tas1-rys1<pys1Cut || tass-ryss<pyssCut || taf-ry<pyallCut  ){
00912         SelectTracks = false;
00913       }
00914       else{
00915         ++niteration;
00916       }
00917       
00918     } while(SelectTracks && niteration < nitMax );      
00919     //======================    finish do loop finder for  xytypecurrent     ====================================================
00920     
00921     //============
00922     
00923     //===========================================================================================================================
00924     
00925     //===========================================================================================================================
00926   }// for xytypecurrent 
00927   //===========================================================================================================================
00928   
00929   if (verbos > 0) {
00930     std::cout << " numberXtracks= " << numberXtracks << " numberYtracks= " << numberYtracks << std::endl;
00931   }
00932   //===========================================================================================================================
00933   //===========================================================================================================================
00934   //===========================================================================================================================
00935 
00936   // case X and Y plane types are available
00937   if(xytype>2) {
00938   //===========================================================================================================================
00939   // match selected X and Y tracks to each other: tgphi=By/Bx->phi=artg(By/Bx); tgtheta=Bx/cosphi=By/sinphi->  ================
00940   //                min of |Bx/cosphi-By/sinphi|                                                               ================
00941 
00942   //  
00943     if (verbos > 0) {
00944       std::cout << " numberXtracks= " << numberXtracks << " numberYtracks= " << numberYtracks << std::endl;
00945     }
00946       if(numberXtracks>0) {
00947         int newxnum[10], newynum[10];// max # tracks = restracks = 10
00948         int nmathed=0;
00949         do {
00950           double dthmin= 999999.; 
00951           int trminx=-1, trminy=-1;
00952           for (int trx=0; trx<numberXtracks; ++trx) {
00953             if (verbos > 0) {
00954               std::cout << "----------- trx= " << trx << " nmathed= " << nmathed << std::endl;
00955             }
00956             for (int tr=0; tr<numberYtracks; ++tr) {
00957               if (verbos > 0) {
00958                 std::cout << "--- tr= " << tr << " nmathed= " << nmathed << std::endl;
00959               }
00960               bool YesY=false;
00961               for (int nmx=0; nmx<nmathed; ++nmx) {
00962                 if(trx==newxnum[nmx]) YesY=true;
00963                 if(YesY) break;
00964                 for (int nm=0; nm<nmathed; ++nm) {
00965                   if(tr==newynum[nm]) YesY=true;
00966                   if(YesY) break;
00967                 }
00968               }
00969               if(!YesY) {
00970                 //--------------------------------------------------------------------  ----    ----    ----    ----    ----    ----
00971                 //double yyyyyy = 999999.;
00972                 //if(Bx[trx] != 0.) yyyyyy = Ay[tr]-(Ax[trx]-xxxvtx)*By[tr]/Bx[trx];
00973                 //double xxxxxx = 999999.;
00974                 //if(By[tr] != 0.) xxxxxx = Ax[trx]-(Ay[tr]-yyyvtx)*Bx[trx]/By[tr];
00975                 //double  dthdif= std::abs(yyyyyy-yyyvtx) + std::abs(xxxxxx-xxxvtx);
00976                 
00977                 double  dthdif= std::abs(AxW[trx]-Ay[tr]) + std::abs(BxW[trx]-By[tr]);
00978                 
00979                 if (verbos > 0) {
00980                   //  std::cout << " yyyyyy= " << yyyyyy << " xxxxxx= " << xxxxxx << " dthdif= " << dthdif << std::endl;
00981                   std::cout << " abs(AxW[trx]-Ay[tr]) = " << std::abs(AxW[trx]-Ay[tr]) << " abs(BxW[trx]-By[tr])= " << std::abs(BxW[trx]-By[tr]) << " dthdif= " << dthdif << std::endl;
00982                 }
00983                 //--------------------------------------------------------------------      ----        ----    ----    ----    ----    ----
00984                 if( dthdif < dthmin ) {
00985                   dthmin = dthdif;
00986                   trminx = trx;
00987                   trminy = tr;
00988                 }//if  dthdif
00989                 //--------------------------------------------------------------------  
00990               }//if !YesY
00991             }//for y
00992           }// for x
00993           ++nmathed;
00994           if(trminx != -1) {
00995             newxnum[nmathed-1] = trminx;
00996           }
00997           else{
00998             newxnum[nmathed-1] = nmathed-1;
00999           }
01000           if (verbos > 0) {
01001             std::cout << " trminx= " << trminx << std::endl;
01002           }
01003           if(nmathed>numberYtracks){
01004             newynum[nmathed-1] = -1;
01005             if (verbos > 0) {
01006               std::cout << "!!!  nmathed= " << nmathed << " > numberYtracks= " << numberYtracks << std::endl;
01007             }
01008           }
01009           else {
01010             if (verbos > 0) {
01011               std::cout << " trminy= " << trminy << std::endl;
01012             }
01013             newynum[nmathed-1] = trminy;
01014           }    
01015         } while(nmathed<numberXtracks && nmathed < restracks);      
01016         
01017         //
01018 //===========================================================================================================================
01019 //
01020     for (int tr=0; tr<nmathed; ++tr) {
01021       int tx=newxnum[tr];
01022       int ty=newynum[tr];
01023       if(ty==-1){
01024         ty=tx;
01025         Ay[ty]=999.;
01026         By[ty]=999.;
01027         Cy[ty]=999.;
01028         My[ty]=-1;
01029       }//if ty
01030       // test:
01031       //  tx=tr;
01032       //ty=tr;
01033       if (verbos > 0) {
01034         if(Mx[tx]>20) {
01035           std::cout << " for track tr= " << tr << " tx= " << tx << " ty= " << ty << std::endl;
01036           std::cout << " Ax= " << Ax[tx]   << " Ay= " << Ay[ty]   << std::endl;
01037           std::cout << " Bx= " << Bx[tx]   << " By= " << By[ty]   << std::endl;
01038           std::cout << " Cx= " << Cx[tx]   << " Cy= " << Cy[ty]   << std::endl;
01039           std::cout << " Mx= " << Mx[tx]   << " My= " << My[ty]   << std::endl;
01040           std::cout << " AxW= " << AxW[tx]   << " AyW= " << AyW[ty]   << std::endl;
01041           std::cout << " BxW= " << BxW[tx]   << " ByW= " << ByW[ty]   << std::endl;
01042           std::cout << " CxW= " << CxW[tx]   << " CyW= " << CyW[ty]   << std::endl;
01043           std::cout << " MxW= " << MxW[tx]   << " MyW= " << MyW[ty]   << std::endl;
01044         }
01045       }
01046       //   rhits.push_back( TrackFP420(c0X,c1X,chisqX,nhitplanesY,c0Y,c1Y,chisqY,nhitplanesY) );
01047       rhits.push_back( TrackFP420(Ax[tx],Bx[tx],Cx[tx],Mx[tx],Ay[ty],By[ty],Cy[ty],My[ty]) );
01048     }//for tr
01049     //============================================================================================================
01050       }//in  numberXtracks >0
01051       //============
01052       
01053   }
01054   // case Y plane types are available only
01055   else if(xytype==1) {
01056     for (int ty=0; ty<numberYtracks; ++ty) {
01057       if (verbos > 0) {
01058         std::cout << " for track ty= " << ty << std::endl;
01059         std::cout << " Ay= " << Ay[ty]   << std::endl;
01060         std::cout << " By= " << By[ty]   << std::endl;
01061         std::cout << " Cy= " << Cy[ty]   << std::endl;
01062         std::cout << " My= " << My[ty]   << std::endl;
01063         std::cout << " AyW= " << AyW[ty]   << std::endl;
01064         std::cout << " ByW= " << ByW[ty]   << std::endl;
01065         std::cout << " CyW= " << CyW[ty]   << std::endl;
01066         std::cout << " MyW= " << MyW[ty]   << std::endl;
01067       }
01068       rhits.push_back( TrackFP420(AyW[ty],ByW[ty],CyW[ty],MyW[ty],Ay[ty],By[ty],Cy[ty],My[ty]) );
01069     }//for ty
01070     //============
01071   }
01072   // case X plane types are available only
01073   else if(xytype==2) {
01074     for (int tx=0; tx<numberXtracks; ++tx) {
01075       if (verbos > 0) {
01076         std::cout << " for track tx= " << tx << std::endl;
01077         std::cout << " Ax= " << Ax[tx]   << std::endl;
01078         std::cout << " Bx= " << Bx[tx]   << std::endl;
01079         std::cout << " Cx= " << Cx[tx]   << std::endl;
01080         std::cout << " Mx= " << Mx[tx]   << std::endl;
01081         std::cout << " AxW= " << AxW[tx]   << std::endl;
01082         std::cout << " BxW= " << BxW[tx]   << std::endl;
01083         std::cout << " CxW= " << CxW[tx]   << std::endl;
01084         std::cout << " MxW= " << MxW[tx]   << std::endl;
01085       }
01086       rhits.push_back( TrackFP420(Ax[tx],Bx[tx],Cx[tx],Mx[tx],AxW[tx],BxW[tx],CxW[tx],MxW[tx]) );
01087     }//for tx
01088     //============
01089   }//xytype
01090 
01091 
01092 
01093 
01095 
01096 
01097 
01098   return rhits;
01099   //============
01100 }
01104 
01105
/data/git/CMSSW_5_3_11_patch5/src/RecoRomanPot/RecoFP420/src/TrackProducerFP420.cc
