--- DarthVader/TrackerLevel2/src/TrkLevel2.cpp	2008/01/22 08:55:07	1.43
+++ DarthVader/TrackerLevel2/src/TrkLevel2.cpp	2008/04/11 13:44:38	1.46
@@ -60,13 +60,20 @@
 	ypu[ip]    = 0;   
 
     };
-
+    
 //     TrkParams::SetTrackingMode();
 //     TrkParams::SetPrecisionFactor();
 //     TrkParams::SetStepMin();
     TrkParams::SetMiniDefault();
     TrkParams::SetPFA();
 
+    int ngf = TrkParams::nGF;
+    for(int i=0; i<ngf; i++){
+	xGF[i] = 0.;
+	yGF[i] = 0.;
+    }
+
+
 };
 //--------------------------------------
 //
@@ -112,6 +119,11 @@
     TrkParams::SetMiniDefault();
     TrkParams::SetPFA();
 
+    int ngf = TrkParams::nGF;
+    for(int i=0; i<ngf; i++){
+	xGF[i] = t.xGF[i];
+	yGF[i] = t.yGF[i];
+    }
 };
 //--------------------------------------
 //
@@ -152,6 +164,11 @@
 	t.ypu[ip]      = ypu[ip];   
 	    
     };
+    int ngf = TrkParams::nGF;
+    for(int i=0; i<ngf; i++){
+	t.xGF[i] = xGF[i];
+	t.yGF[i] = yGF[i];
+    }
 
     
 };
@@ -165,9 +182,15 @@
  * @param t pointer to an object of the class Trajectory, 
  * which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
  * @return error flag.
+ *
+ * >>> OBSOLETE !!! use TrkTrack::DoTrack2(Trajectory* t) instead
+ *
  */
 int TrkTrack::DoTrack(Trajectory* t){
 
+    cout << " int TrkTrack::DoTrack(Trajectory* t) --->> OBSOLETE !!! "<<endl;
+    cout << " use int TrkTrack::DoTrack2(Trajectory* t)"<<endl;
+
     double *dxout = new double[t->npoint];
     double *dyout = new double[t->npoint];
     double *dzin = new double[t->npoint];
@@ -186,13 +209,13 @@
     dotrack_(&(t->npoint),dzin,dxout,dyout,dal,&ifail);
     
     for (int i=0; i<t->npoint; i++){
-	t->x[i] = (float)*dxout++;
-	t->y[i] = (float)*dyout++;
+	t->x[i] = (float)*(dxout+i);
+	t->y[i] = (float)*(dyout+i);
     }
 
-//    delete [] dxout;
-//    delete [] dyout;
-//    delete [] dzin;
+    delete [] dxout;
+    delete [] dyout;
+    delete [] dzin;
 
     return ifail;
 };
@@ -230,16 +253,19 @@
     dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
     
     for (int i=0; i<t->npoint; i++){
-	t->x[i]   = (float)*dxout++;
-	t->y[i]   = (float)*dyout++;
-	t->thx[i] = (float)*dthxout++;
-	t->thy[i] = (float)*dthyout++;
-	t->tl[i]  = (float)*dtlout++;
+	t->x[i]   = (float)*(dxout+i);
+	t->y[i]   = (float)*(dyout+i);
+	t->thx[i] = (float)*(dthxout+i);
+	t->thy[i] = (float)*(dthyout+i);
+	t->tl[i]  = (float)*(dtlout+i);
     }
 
-//    delete [] dxout;
-//    delete [] dyout;
-//    delete [] dzin;
+    delete [] dxout;
+    delete [] dyout;
+    delete [] dzin;
+    delete [] dthxout;
+    delete [] dthyout;
+    delete [] dtlout;
 
     return ifail;
 };
@@ -299,7 +325,8 @@
     return dedx;
 };
 /**
- * Returns 1 if the cluster on a tracker view includes bad strips.
+ * Returns 1 if the cluster on a tracker view includes bad strips 
+ * (at least one bad strip among the four strip used by p.f.a.)
  * @param ip plane (0-5)
  * @param iv view (0=x 1=y)
  */
@@ -651,22 +678,64 @@
 	track.zm[i]=zm[i];
 	
 //      --- temporaneo ----------------------------
-//      andrebbe inserita la dimensione del sensore
-	float segment = 100.;
-	track.xm_a[i]=xm[i];
-	track.xm_b[i]=xm[i];
-	track.ym_a[i]=ym[i];
-	track.ym_b[i]=ym[i];
-	if(       XGood(i) && !YGood(i) ){
-	    track.ym_a[i] = track.ym_a[i]+segment;
-	    track.ym_b[i] = track.ym_b[i]-segment;
-	}else if( !XGood(i) && YGood(i)){
-	    track.xm_a[i] = track.xm_a[i]+segment;
-	    track.xm_b[i] = track.xm_b[i]-segment;
-	}
+// 	float segment = 100.;
+// 	track.xm_a[i]=xm[i];
+// 	track.xm_b[i]=xm[i];
+// 	track.ym_a[i]=ym[i];
+// 	track.ym_b[i]=ym[i];
+// 	if(       XGood(i) && !YGood(i) ){
+// 	    track.ym_a[i] = track.ym_a[i]+segment;
+// 	    track.ym_b[i] = track.ym_b[i]-segment;
+// 	}else if( !XGood(i) && YGood(i)){
+// 	    track.xm_a[i] = track.xm_a[i]+segment;
+// 	    track.xm_b[i] = track.xm_b[i]-segment;
+// 	}
 //      --- temporaneo ----------------------------
+
+	if( XGood(i) || YGood(i) ){
+	    double segment = 2.;//cm
+	    // NB: i parametri di allineamento hanno una notazione particolare!!!
+	    // sensor = 0 (hybrid side), 1 
+	    // ladder = 0-2 (increasing x)
+	    // plane  = 0-5 (from bottom to top!!!) 
+	    int is = (int)GetSensor(i); if(i==5)is=1-is;
+	    int ip = 5-i;
+	    int il = (int)GetLadder(i);
+	    
+	    double omega   = 0.; 
+	    double beta    = 0.;
+	    double gamma   = 0.;
+	    if( 
+		(is < 0 || is > 1 || ip < 0 || ip > 5 || il < 0 || il > 2) &&
+		true){
+		// se il piano risulta colpito, ladder e sensore devono essere 
+		// assegnati correttamente
+		cout << " void TrkTrack::FillMiniStruct(cMini2track&) --- WARNING --- sensor not defined, cannot read alignment parameters "<<endl;
+		cout << " is ip il = "<<is<<" "<<ip<<" "<<il<<endl;
+	    }else{
+		omega   = alignparameters_.omega[is][il][ip];
+		beta    = alignparameters_.beta[is][il][ip];
+		gamma   = alignparameters_.gamma[is][il][ip];
+	    }
+	    
+	    if(       XGood(i) && !YGood(i) ){
+		track.xm_a[i] = xm[i] - omega * segment;
+		track.ym_a[i] = ym[i] + segment;
+//	    track.zm_a[i] = zm[i] + beta * segment;//not used yet
+		track.xm_b[i] = xm[i] + omega * segment;
+		track.ym_b[i] = ym[i] - segment;
+//	    track.zm_b[i] = zm[i] - beta * segment;//not used yet
+	    }else if( !XGood(i) && YGood(i) ){
+		track.xm_a[i] = xm[i] + segment;
+		track.ym_a[i] = ym[i] + omega * segment;
+//	    track.zm_a[i] = zm[i] - gamma * segment;//not used yet
+		track.xm_b[i] = xm[i] - segment;
+		track.ym_b[i] = ym[i] - omega * segment;
+//	    track.zm_b[i] = zm[i] + gamma * segment;//not used yet
+	    }
+	}
 	
-	track.resx[i]=resx[i];
+       	track.resx[i]=resx[i];
 	track.resy[i]=resy[i];
 	track.tailx[i]=tailx[i];
 	track.taily[i]=taily[i];
@@ -890,20 +959,69 @@
 }
 
 /**
- * Returns 1 if the track is inside the magnet cavity
- * Set the minimum number of steps for tracking precision 
+ * Returns true if the track is inside the magnet cavity.
+ * @param toll Tolerance around the nominal volume (toll>0 define an inner fiducial volume)
  */
-Bool_t TrkTrack::IsInsideCavity(){
-    float xmagntop, ymagntop, xmagnbottom, ymagnbottom;
-    xmagntop = xv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*axv[0]/180.);
-    ymagntop = yv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*ayv[0]/180.);
-    xmagnbottom = xv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*axv[5]/180.);
-    ymagnbottom = yv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*ayv[5]/180.);
-    if( xmagntop>XMAGNLOW && xmagntop<XMAGNHIGH &&
-	ymagntop>YMAGNLOW && ymagntop<YMAGNHIGH &&
-	xmagnbottom>XMAGNLOW && xmagnbottom<XMAGNHIGH &&
-	ymagnbottom>YMAGNLOW && ymagnbottom<YMAGNHIGH ) return(true);
-    else return(false);
+Bool_t TrkTrack::IsInsideCavity(float toll){
+
+//     float xmagntop, ymagntop, xmagnbottom, ymagnbottom;
+//     xmagntop = xv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*axv[0]/180.);
+//     ymagntop = yv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*ayv[0]/180.);
+//     xmagnbottom = xv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*axv[5]/180.);
+//     ymagnbottom = yv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*ayv[5]/180.);
+//     if( xmagntop>XMAGNLOW && xmagntop<XMAGNHIGH &&
+// 	ymagntop>YMAGNLOW && ymagntop<YMAGNHIGH &&
+// 	xmagnbottom>XMAGNLOW && xmagnbottom<XMAGNHIGH &&
+// 	ymagnbottom>YMAGNLOW && ymagnbottom<YMAGNHIGH ) return(true);
+//     else return(false);
+
+    int ngf = TrkParams::nGF;
+    for(int i=0; i<ngf; i++){
+	//
+//	cout << endl << TrkParams::GF_element[i];
+	if( 
+	    TrkParams::GF_element[i].CompareTo("CUF") &&
+	    TrkParams::GF_element[i].CompareTo("T2")  &&
+	    TrkParams::GF_element[i].CompareTo("T3")  &&
+	    TrkParams::GF_element[i].CompareTo("T4")  &&
+	    TrkParams::GF_element[i].CompareTo("T5")  &&
+	    TrkParams::GF_element[i].CompareTo("CLF") &&
+	    true)continue;
+	// apply condition only within the cavity
+//	cout << " -- "<<xGF[i]<<" "<<yGF[i];
+	if( 
+	    xGF[i] <= TrkParams::xGF_min[i] + toll ||
+	    xGF[i] >= TrkParams::xGF_max[i] - toll ||
+	    yGF[i] <= TrkParams::yGF_min[i] + toll ||
+	    yGF[i] >= TrkParams::yGF_max[i] - toll ||
+	    false){
+	    
+	    return false;
+	}
+    }
+    return true;
+
+
+}
+/**
+ * Returns true if the track is inside the nominal acceptance, which is defined 
+ * by the intersection among magnet cavity, silicon-plane sensitive area and 
+ * ToF sensitive area (nominal values from the official document used to 
+ * calculate the geometrical factor)
+ */
+Bool_t TrkTrack::IsInsideAcceptance(){
+
+    int ngf = TrkParams::nGF;
+    for(int i=0; i<ngf; i++){
+	if( 
+	    xGF[i] <= TrkParams::xGF_min[i] ||
+	    xGF[i] >= TrkParams::xGF_max[i] ||
+	    yGF[i] <= TrkParams::yGF_min[i] ||
+	    yGF[i] >= TrkParams::yGF_max[i] ||
+	    false)return false;
+    }
+    return true;
+
 }
 /**
  * Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track.
@@ -923,7 +1041,7 @@
 };
 
 /**
- * Method to retrieve the ladder (0-4, increasing x) traversed by the track on this plane.
+ * Method to retrieve the ladder (0-2, increasing x) traversed by the track on this plane.
  * If no ladder is traversed (dead area) the metod retuns -1.
  * @param ip Tracker plane (0-5)
  */
@@ -946,27 +1064,35 @@
 /**
  * \brief Method to include a x-cluster to the track.
  * @param ip Tracker plane (0-5)
- * @param clid Cluster ID (0,1,...)
- * @param is Sensor (0-1, increasing y) 
+ * @param clid Cluster ID (0 = no-cluster, 1,2,... otherwise )
+ * @param il Ladder (0-2, increasing x, -1 if no sensitive area is hit) 
+ * @param is Sensor (0-1, increasing y, -1 if no sensitive area is hit) 
+ * @param bad True if the cluster contains bad strips  
  * @see Fit(double pfixed, int& fail, int iprint, int froml1)
  */
-void TrkTrack::SetXGood(int ip, int clid, int is){
-    int il=0;       //ladder (temporary)
-    bool bad=false; //ladder (temporary)
-    xgood[ip]=il*100000000+is*10000000+clid;
+void TrkTrack::SetXGood(int ip, int clid, int il, int is, bool bad){
+//    int il=0;       //ladder (temporary)
+//    bool bad=false; //ladder (temporary)
+    if(ip<0||ip>5||clid<0||il<-1||il>2||is<-1||is>1)
+	cout << " void TrkTrack::SetXGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl;
+    xgood[ip]=(il+1)*100000000+(is+1)*10000000+clid;
     if(bad)xgood[ip]=-xgood[ip];
 };
 /**
  * \brief Method to include a y-cluster to the track.
  * @param ip Tracker plane (0-5)
- * @param clid Cluster ID (0,1,...)
- * @param is Sensor (0-1)
+ * @param clid Cluster ID (0 = no-cluster, 1,2,... otherwise )
+ * @param il Ladder (0-2, increasing x, -1 if no sensitive area is hit) 
+ * @param is Sensor (0-1, increasing y, -1 if no sensitive area is hit) 
+ * @param bad True if the cluster contains bad strips  
  * @see Fit(double pfixed, int& fail, int iprint, int froml1)
  */
-void TrkTrack::SetYGood(int ip, int clid, int is){
-    int il=0;       //ladder (temporary)
-    bool bad=false; //ladder (temporary)
-    ygood[ip]=il*100000000+is*10000000+clid;
+void TrkTrack::SetYGood(int ip, int clid, int il, int is, bool bad){
+//    int il=0;       //ladder (temporary)
+//    bool bad=false; //ladder (temporary)
+    if(ip<0||ip>5||clid<0||il<-1||il>2||is<-1||is>1)
+	cout << " void TrkTrack::SetYGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl;
+    ygood[ip]=(il+1)*100000000+(is+1)*10000000+clid;
     if(bad)ygood[ip]=-ygood[ip];
 };
 
@@ -1219,6 +1345,11 @@
 	dedx_y[ip] = 0;
 
     };
+    int ngf = TrkParams::nGF;
+    for(int i=0; i<ngf; i++){
+	xGF[i] = 0.;
+	yGF[i] = 0.;
+    }
 //     if(clx)clx->Clear();
 //     if(cly)cly->Clear();
 //    clx.Clear();
@@ -1245,11 +1376,13 @@
 //--------------------------------------
 TrkSinglet::TrkSinglet(){
 //    cout << "TrkSinglet::TrkSinglet() " << GetUniqueID()<<endl;
-    plane    = 0;
-    coord[0] = 0; 
-    coord[1] = 0;
-    sgnl     = 0;
+//     plane    = 0;
+//     coord[0] = 0; 
+//     coord[1] = 0;
+//     sgnl     = 0;
+//     multmax  = 0;
 //    cls      = 0;
+    Clear();
 };
 //--------------------------------------
 //
@@ -1261,6 +1394,7 @@
     coord[0] = s.coord[0]; 
     coord[1] = s.coord[1];
     sgnl     = s.sgnl;
+    multmax  = s.multmax;
 //	cls      = 0;//<<<<pointer
 //    cls      = TRef(s.cls);
 };
@@ -1271,9 +1405,11 @@
 void TrkSinglet::Dump(){
     int i=0;
     cout << endl << "========== Singlet " ;
-    cout << endl << "plane    : " << plane;
-    cout << endl << "coord[2] : "; while( i<2 && cout << coord[i] << " ") i++;
-    cout << endl << "sgnl     : " << sgnl;
+    cout << endl << "plane        : " << plane;
+    cout << endl << "coord[2]     : "; while( i<2 && cout << coord[i] << " ") i++;
+    cout << endl << "sgnl         : " << sgnl;
+    cout << endl << "max.strip    : ";
+    cout << endl << "multiplicity : ";
 }
 //--------------------------------------
 //
@@ -1286,6 +1422,7 @@
     coord[0]=-999;
     coord[1]=-999;
     sgnl=0;
+    multmax  = 0;
     
 }
 //--------------------------------------
@@ -1396,7 +1533,7 @@
  * forced (see TrkLevel2::GetVKMask(int,int)) or 
  * for this event only (TrkLevel2::GetVKFlag(int,int)). 
  * @param iv Tracker view (0-11)
- * @param ivk Viking-chip number (0-23)
+ * @param ivk Viking-chip number (0-23) 
  */
 Bool_t TrkLevel2::IsMaskedVK(int iv, int ivk){
     return !(GetVKMask(iv,ivk)&&GetVKFlag(iv,ivk) ); 
@@ -1492,10 +1629,22 @@
 	    //-----------------------------------------------------
 	    //-----------------------------------------------------
 	};
+	// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+	// evaluated coordinates (to define GF)
+	// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+	int    ngf = TrkParams::nGF;
+	float *zgf = TrkParams::zGF;
+	Trajectory tgf = Trajectory(ngf,zgf);
+	tgf.DoTrack2(t_track->al);//<<<< integrate the trajectory
+	for(int ip=0; ip<ngf; ip++){
+	    t_track->xGF[ip] = tgf.x[ip];
+	    t_track->yGF[ip] = tgf.y[ip];
+	}
+	
 // 	if(t_track->IsSaturated())t_track->Dump();
 	new(t[i]) TrkTrack(*t_track);
 	t_track->Clear();
-    };
+    };//end loop over track
 
 //  ----------------
 //  *** SINGLETS ***
@@ -1507,6 +1656,8 @@
 	t_singlet->coord[0] = l2->xs[i][0];
 	t_singlet->coord[1] = l2->xs[i][1];
 	t_singlet->sgnl     = l2->signlxs[i];
+	t_singlet->multmax = l2->multmaxsx[i]; 
+	if(l2->sxbad[i]>0) t_singlet->multmax = -1*t_singlet->multmax;
 	//-----------------------------------------------------
 //	if(l1) t_singlet->cls      = l1->GetCluster(l2->clsx[i]-1);
 	//-----------------------------------------------------
@@ -1520,12 +1671,16 @@
 	t_singlet->coord[0] = l2->ys[i][0];
 	t_singlet->coord[1] = l2->ys[i][1];
 	t_singlet->sgnl     = l2->signlys[i];
+	t_singlet->multmax  = l2->multmaxsy[i]; 
+	if(l2->sybad[i]>0) t_singlet->multmax = -1*t_singlet->multmax;
 	//-----------------------------------------------------
 //	if(l1) t_singlet->cls      = l1->GetCluster(l2->clsy[i]-1);
 	//-----------------------------------------------------
 	new(sy[i]) TrkSinglet(*t_singlet);
 	t_singlet->Clear();
     };
+
+
 	
     delete t_track;
     delete t_singlet;
@@ -2040,13 +2195,28 @@
  */
 int Trajectory::DoTrack2(float* al){
 
-    double *dxout   = new double[npoint];
-    double *dyout   = new double[npoint];
-    double *dthxout = new double[npoint];
-    double *dthyout = new double[npoint];
-    double *dtlout  = new double[npoint];
-    double *dzin    = new double[npoint];
-    double dal[5];
+//      double *dxout   = new double[npoint];
+//      double *dyout   = new double[npoint];
+//      double *dthxout = new double[npoint];
+//      double *dthyout = new double[npoint];
+//      double *dtlout  = new double[npoint];
+//      double *dzin    = new double[npoint];
+    
+     double *dxout;
+     double *dyout;
+     double *dthxout;
+     double *dthyout;
+     double *dtlout;
+     double *dzin;
+
+     dxout   = (double*) malloc(npoint*sizeof(double));
+     dyout   = (double*) malloc(npoint*sizeof(double));
+     dthxout = (double*) malloc(npoint*sizeof(double));
+     dthyout = (double*) malloc(npoint*sizeof(double));
+     dtlout  = (double*) malloc(npoint*sizeof(double));
+     dzin    = (double*) malloc(npoint*sizeof(double));
+
+     double dal[5];
 
     int ifail = 0;
 
@@ -2061,12 +2231,27 @@
     dotrack2_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
     
     for (int i=0; i<npoint; i++){
-	x[i]   = (float)*dxout++;
-	y[i]   = (float)*dyout++;
-	thx[i] = (float)*dthxout++;
-	thy[i] = (float)*dthyout++;
-	tl[i]  = (float)*dtlout++;
-    }
+	x[i]   = (float)*(dxout+i);
+	y[i]   = (float)*(dyout+i);
+	thx[i] = (float)*(dthxout+i);
+	thy[i] = (float)*(dthyout+i);
+	tl[i]  = (float)*(dtlout+i);
+    }
+
+    if(dxout)  free( dxout );
+    if(dyout)  free( dyout );
+    if(dthxout)free( dthxout );
+    if(dthyout)free( dthyout );
+    if(dtlout) free( dtlout );
+    if(dzin)   free( dzin );
+
+//      delete [] dxout;
+//      delete [] dyout;
+//      delete [] dthxout;
+//      delete [] dthyout;
+//      delete [] dtlout;
+//      delete [] dzin;
+
 
     return ifail;
 };