FairRoot/PandaRoot
PndTrkGlpkFits.cxx
Go to the documentation of this file.
1 #include "PndTrkGlpkFits.h"
2 #include "glpk.h"
3 #include <cmath>
4 #include <iostream>
5 // Root includes
6 #include "TROOT.h"
7 
8 using namespace std;
9 
10 #define PI 3.141592654
11 
12 //----------begin of function PndTrkGlpkFits::FitHelixCylinder
13 
14 Short_t PndTrkGlpkFits::FitHelixCylinder(
15  Short_t nHitsinTrack,
16  Double_t *Xconformal,
17  Double_t *Yconformal,
18  Double_t *DriftRadiusconformal,
19  Double_t *ErrorDriftRadiusconformal,
20  Double_t rotationangle,
21  Double_t trajectory_vertex[2],
22  Short_t NMAX,
23  Double_t *emme,
24  Double_t *qu,
25  Double_t *pAlfa,
26  Double_t *pBeta,
27  Double_t *pGamma,
28  bool *Type,
29  int istampa,
30  int IVOLTE
31  )
32 {
33 
34 
35  // definition of variables for the glpsol solver
36  // ROWS (for read_rows function)
37  //
38 
39  Short_t NpointsInFit = nHitsinTrack-NMAX <0 ? nHitsinTrack : NMAX;
40  bool mvdhit[NpointsInFit];
41 
42  Double_t
43  A,
44  alfetta,
45  angle,
46  Delta[NpointsInFit],
47  M = 1.,
48  m_result,
49  offsety,
50  Oxx[NpointsInFit],
51  Oyy[NpointsInFit],
52  q_result;
53 
54  Short_t i, j, ii, iii, nSttHits, nMvdHits;
55  Short_t Status;
56 
57  float m1_result,m2_result, q1_result,q2_result, A1_result, A2_result;
58 
59 // --
60 
61 //-------------- stampaggi
62 if(istampa>=3){
63  cout<<"from FitHelixCylinder,prima di rotazione, Evento "<<IVOLTE<<", nHitsinTrack = "<<nHitsinTrack
64  <<"\nfrom FitHelixCylinder, nPointsinFit = "<<NpointsInFit<<endl;
65  for(i=0 ; i< NpointsInFit ; i++) {
66  cout<<" Xconformal["<<i<<"] = "<<Xconformal[ i ]<<
67  "; Yconformal["<<i<<"] = "<<Yconformal[ i ]<<", drift radius conformal "<<
68  DriftRadiusconformal[i]<<endl<<"\tErrordiriftradiusconformal = "
69  <<ErrorDriftRadiusconformal[i]<<endl;
70  }
71 }
72 //------------ end stampaggi
73 
74 
75 
76  if( nHitsinTrack < 2) {
77  *Type = false ;
78  return -1;
79  }
80 // use the trick of increasing the rotation angle by 10 degrees in order to obtain always a positive m
81  rotationangle -= PI/18.;
82 
83  Double_t cose = cos(rotationangle), sine = sin(rotationangle);
84 
85 
86 
87  nSttHits = nMvdHits = 0;
88  for(i=0;i<NpointsInFit; i++){
89  Oxx[i] = Xconformal[ i ] *cose +
90  Yconformal[ i ]*sine;
91  Oyy[i] = -Xconformal[ i ] *sine +
92  Yconformal[ i ]*cose;
93  Delta[i] = 3.*ErrorDriftRadiusconformal[ i ]; // 3 times the Drift Radius
94 
95  if( DriftRadiusconformal[ i ]<0. )
96  {
97  mvdhit[i]=true;
98  nMvdHits++;
99  } else {
100  mvdhit[i]=false;
101  nSttHits++;
102  }
103  }
104 
105 //-------------- stampaggi
106 if(istampa>=3){
107  cout<<"from FitHelixCylinder, dopo rotazione, Evento "<<IVOLTE
108  <<", nHitsinTrack = "<<nHitsinTrack
109  <<"\nfrom FitHelixCylinder, nPointsinFit = "<<NpointsInFit<<endl;
110  for(i=0 ; i< NpointsInFit ; i++) {
111  cout<<" Ox["<<i<<"] = "<<Oxx[ i ]<<
112  "; Oy["<<i<<"] = "<<Oyy[ i ]<<", Delta "<<
113  Delta[i]<<endl;
114  }
115 }
116 //------------ end stampaggi
117 //--------- calculation of # structural variables (see Gianluigi's logbook pag. 236) etc.etc.
118 
119  int NStructVar = 4 + 1 + nMvdHits * 2 + nSttHits *4 ;
120 // m1,m2,q1,q2 DUM lam & sigma lamp & lamm & sigmap & sigmam
121 
122 
123  int nRows = 1 + nMvdHits * 4 + nSttHits * 9;
124 // OBJECT A,B,C,D Ap,Bp,Cp,Dp,Am,Bm,Cm,Dm,LM
125  int NStructRowsMax = 8*NpointsInFit ; // maximum number of ROWS in which a
126  // structural variable (for instance M ) can be found
127  int NRowsInWhichStructVarArePresent[NStructVar];
128  int nRanges = nSttHits;
129  int nBounds=NpointsInFit+nSttHits+1;
130 
131 
132 
133 //---- creating the various service arrays
134 
135  char OBJECTname[8];
136  sprintf(OBJECTname,"OBJECT");
137 
138  int typeRows[nRows];
139  char * nameRows[nRows];
140  char auxnameRows[nRows][5];
141 
142  double final_values[NStructVar];
143  char *StructVarName[NStructVar];
144  char auxStructVarName[NStructVar][8];
145  char *NameRowsInWhichStructVarArePresent[NStructVar*NStructRowsMax];
146  char aux[NStructVar*NStructRowsMax][5];
147  double Coefficients[NStructVar*NStructRowsMax];
148  //--------for RHS information
149  double ValueB[nRows-1]; // -1 because OBJECT dowsn't have a RHS boundary.
150  //--------for RANGES information
151  double ValueRanges[nRanges];
152  char* NameRanges[nRanges];
153  char auxNameRanges[nRanges][20];
154 
155  //--------for BOUNDS information
156  double BoundValue[nBounds];
157  char *BoundStructVarName[nBounds];
158  char auxBoundStructVarName[nBounds][8];
159  char *TypeofBound[nBounds];
160  char auxTypeofBound[nBounds][20];
161  //--------end BOUNDS information
162 
163 //---------------------------------------------------
164 
165 //--- calculate array NRowsInWhichStructVarArePresent
166  NRowsInWhichStructVarArePresent[0] = // this is for m1
167  NRowsInWhichStructVarArePresent[1] = // this is for m2
168  NRowsInWhichStructVarArePresent[2] = // this is for q1
169  NRowsInWhichStructVarArePresent[3] = nMvdHits*2 + nSttHits *4; // this is for q2
170  //--- the following is for the lam* (Mvd hits) or lamp* (Stt hits) structural variables
171  ii=0;
172  for(i=0; i< NpointsInFit ; i++) {
173  if( mvdhit[i]){
174  NRowsInWhichStructVarArePresent[4+ii]= 4;
175  } else {
176  NRowsInWhichStructVarArePresent[4+ii]= 5;
177  }
178  ii++;
179  }
180  //--- the following is for the lamm* (Mvd Hits, if any) structural variables
181  for(i=0; i< NpointsInFit ; i++) {
182  if( !mvdhit[i]){
183  NRowsInWhichStructVarArePresent[4+ii]= 5;
184  ii++;
185  }
186  }
187 
188  //--- the following is for the sigma structural variables
189  for(i=0 ; i< nMvdHits+2*nSttHits ; i++) {
190  NRowsInWhichStructVarArePresent[4+ii+i]= 5;
191  }
192 //--- the following is for the DUM structural variable
193 
194  NRowsInWhichStructVarArePresent[4+ii+nMvdHits+2*nSttHits]= nMvdHits*4 + nSttHits*8;
195 
196 
197 //----------------- write the ROWS section
198 // sprintf(&(auxnameRows[0][0]),"OBJECT");
199 // nameRows[0]=&auxnameRows[0][0];
200  nameRows[0]=OBJECTname;
201  typeRows[0]=GLP_FR;
202  for(i=0 , ii=0 ; i< NpointsInFit ; i++) {
203 
204  if( mvdhit[i]){
205  typeRows[1+ii]=GLP_UP;typeRows[2+ii]=GLP_UP;typeRows[3+ii]=GLP_UP;typeRows[4+ii]=GLP_UP;
206  typeRows[5+ii]=GLP_LO;
207 
208  sprintf(&(auxnameRows[1+ii][0]),"A%d",i); nameRows[1+ii]=&auxnameRows[1+ii][0];
209  sprintf(&(auxnameRows[2+ii][0]),"B%d",i); nameRows[2+ii]=&auxnameRows[2+ii][0];
210  sprintf(&(auxnameRows[3+ii][0]),"C%d",i); nameRows[3+ii]=&auxnameRows[3+ii][0];
211  sprintf(&(auxnameRows[4+ii][0]),"D%d",i); nameRows[4+ii]=&auxnameRows[4+ii][0];
212  ii+=4;
213 
214  } else {
215  typeRows[1+ii]=GLP_UP;typeRows[2+ii]=GLP_UP;typeRows[3+ii]=GLP_UP;typeRows[4+ii]=GLP_UP;
216  typeRows[5+ii]=GLP_UP;typeRows[6+ii]=GLP_UP;typeRows[7+ii]=GLP_UP;typeRows[8+ii]=GLP_UP;
217  typeRows[9+ii]=GLP_LO;
218 
219  sprintf(&(auxnameRows[1+ii][0]),"Ap%d",i); nameRows[1+ii]=&auxnameRows[1+ii][0];
220  sprintf(&(auxnameRows[2+ii][0]),"Bp%d",i); nameRows[2+ii]=&auxnameRows[2+ii][0];
221  sprintf(&(auxnameRows[3+ii][0]),"Cp%d",i); nameRows[3+ii]=&auxnameRows[3+ii][0];
222  sprintf(&(auxnameRows[4+ii][0]),"Dp%d",i); nameRows[4+ii]=&auxnameRows[4+ii][0];
223  sprintf(&(auxnameRows[5+ii][0]),"Am%d",i); nameRows[5+ii]=&auxnameRows[5+ii][0];
224  sprintf(&(auxnameRows[6+ii][0]),"Bm%d",i); nameRows[6+ii]=&auxnameRows[6+ii][0];
225  sprintf(&(auxnameRows[7+ii][0]),"Cm%d",i); nameRows[7+ii]=&auxnameRows[7+ii][0];
226  sprintf(&(auxnameRows[8+ii][0]),"Dm%d",i); nameRows[8+ii]=&auxnameRows[8+ii][0];
227  sprintf(&(auxnameRows[9+ii][0]),"LM%d",i); nameRows[9+ii]=&auxnameRows[9+ii][0];
228  ii+=9;
229  }
230  }
231 
232 
233 
234 
235 //----------------- write the COLUMNS section
236 
237 // fprintf(FMCS,"COLUMNS\n"); //--------stampaggi
238 
239 // Column variable m1
240 
241  ii=0;
242  for(i=0 ; i< NpointsInFit ; i++) {
243 
244  if( mvdhit[i]){
245 //---stampaggi
246 // fprintf(FMCS," m1 A%d %g\n m1 B%d %g\n",i,Oxx[i],i,-Oxx[i]);
247 //-----stampaggi, fine
248 
249  Coefficients[ii]= Oxx[i];
250  Coefficients[ii+1]= -Oxx[i];
251  ii +=2;
252  } else {
253 //---stampaggi
254 // fprintf(FMCS," m1 Ap%d %g Am%d %g\n m1 Bp%d %g Bm%d %g\n",
255 // i,Oxx[i],i,Oxx[i],i,-Oxx[i],i,-Oxx[i]);
256 //-----stampaggi, fine
257 
258  Coefficients[ii]= Oxx[i];
259  Coefficients[ii+1]= Oxx[i];
260  Coefficients[ii+2]= -Oxx[i];
261  Coefficients[ii+3]= -Oxx[i];
262  ii += 4;
263  }
264  }
265 
266 
267 
268 // Column variable m2
269  for(i=0, ii=0; i< NpointsInFit ; i++) {
270  if( mvdhit[i]){
271 // fprintf(FMCS," m2 A%d %g\n m2 B%d %g\n",//---stampaggi
272 // i,-Oxx[i],i,Oxx[i]);//---stampaggi
273  Coefficients[NStructRowsMax+ii]= -Oxx[i];
274  Coefficients[NStructRowsMax+ii+1]= Oxx[i];
275  ii += 2;
276  } else {
277 // fprintf(FMCS," m2 Ap%d %g Am%d %g\n m2 Bp%d %g Bm%d %g\n",//---stampaggi
278 // i,-Oxx[i],i,-Oxx[i],i,Oxx[i],i,Oxx[i]);//---stampaggi
279  Coefficients[NStructRowsMax+ii]= -Oxx[i];
280  Coefficients[NStructRowsMax+ii+1]= -Oxx[i];
281  Coefficients[NStructRowsMax+ii+2]= Oxx[i];
282  Coefficients[NStructRowsMax+ii+3]= Oxx[i];
283  ii += 4;
284  }
285  }
286 
287 // Column variable q1
288  for(i=0, ii=0 ; i< NpointsInFit ; i++) {
289  if( mvdhit[i]){
290 // fprintf(FMCS," q1 A%d 1.\n q1 B%d -1.\n",//---stampaggi
291 // i,i);//---stampaggi
292  Coefficients[2*NStructRowsMax+ii]= 1.;
293  Coefficients[2*NStructRowsMax+ii+1]= -1.;
294  ii +=2;
295  } else {
296 // fprintf(FMCS," q1 Ap%d 1. Am%d 1.\n q1 Bp%d -1. Bm%d -1.\n",//---stampaggi
297 // i,i,i,i);//---stampaggi
298  Coefficients[2*NStructRowsMax+ii]= 1.;
299  Coefficients[2*NStructRowsMax+ii+1]= 1.;
300  Coefficients[2*NStructRowsMax+ii+2]= -1.;
301  Coefficients[2*NStructRowsMax+ii+3]= -1.;
302  ii += 4;
303  }
304  }
305 
306 // Column variable q2
307  for(i=0, ii=0 ; i< NpointsInFit ; i++) {
308  if( mvdhit[i]){
309 // fprintf(FMCS," q2 A%d -1.\n q2 B%d 1.\n",//---stampaggi
310 // i,i);//---stampaggi
311  Coefficients[3*NStructRowsMax+ii]= -1.;
312  Coefficients[3*NStructRowsMax+ii+1]= 1.;
313  ii += 2;
314  } else {
315 // fprintf(FMCS," q2 Ap%d -1. Am%d -1.\n q2 Bp%d 1. Bm%d 1.\n",//---stampaggi
316 // i,i,i,i);//---stampaggi
317  Coefficients[3*NStructRowsMax+ii]= -1.;
318  Coefficients[3*NStructRowsMax+ii+1]= -1.;
319  Coefficients[3*NStructRowsMax+ii+2]= 1.;
320  Coefficients[3*NStructRowsMax+ii+3]= 1.;
321  ii += 4;
322  }
323 
324  }
325 
326 // Column variable lambdap(i)
327  for(i=0 ; i< NpointsInFit ; i++) {
328  ii=(4+i)*NStructRowsMax;
329  Coefficients[ii]= -M;
330  Coefficients[ii+1]= -M;
331  Coefficients[ii+2]= -M;
332  Coefficients[ii+3]= M;
333 // if( mvdhit[i]){
334 // fprintf(FMCS," lam%d A%d %g B%d %g\n lam%d C%d %g D%d %g\n",//---stampaggi
335 // i,i,-M,i,-M, i , i,-M, i, M);//---stampaggi
336 // } else {
337 // fprintf(FMCS," lamp%d Ap%d %g Bp%d %g\n lamp%d Cp%d %g Dp%d %g\n lamp%d LM%d 1.\n",//---stampaggi
338 // i,i,-M,i,-M, i , i,-M, i, M, i,i);//---stampaggi
339 // }
340 
341  if(! mvdhit[i]) Coefficients[ii+4]= 1.;
342 
343  }
344 
345 
346 
347 // Column variable lambdam(i)
348  for(i=0 ; i< NpointsInFit ; i++) {
349  if( mvdhit[i]) continue;
350  ii+= NStructRowsMax;
351 // fprintf(FMCS," lamm%d Am%d %g Bm%d %g\n lamm%d Cm%d %g Dm%d %g\n lamm%d LM%d 1.\n",//---stampaggi
352 // i,i,-M,i,-M, i,i , -M, i, M, i,i);//---stampaggi
353  Coefficients[ii]= -M;
354  Coefficients[ii+1]= -M;
355  Coefficients[ii+2]= -M;
356  Coefficients[ii+3]= M;
357  Coefficients[ii+4]= 1.;
358  }
359 // Column variable SIGp(i)
360  for(i=0; i< NpointsInFit ; i++) {
361  ii+= NStructRowsMax;
362 
363 /*
364  if( mvdhit[i]){
365  fprintf(FMCS," SIG%d OBJECT %g A%d -1.\n SIG%d B%d -1. C%d 1.\n SIG%d D%d -1.\n",//---stampaggi
366  i,1./Delta[i],i,i,i,i,i,i);//---stampaggi
367 
368  } else {
369  fprintf(FMCS," SIGp%d OBJECT %g Ap%d -1.\n SIGp%d Bp%d -1. Cp%d 1.\n SIGp%d Dp%d -1.\n",//---stampaggi
370  i,1./Delta[i],i,i,i,i,i,i);//---stampaggi
371  }
372 */
373 
374  Coefficients[ii]= 1./Delta[i];
375  Coefficients[ii+1]= -1.;
376  Coefficients[ii+2]= -1.;
377  Coefficients[ii+3]= 1.;
378  Coefficients[ii+4]= -1.;
379 
380 
381  }
382 // Column variable SIGm(i)
383  for(i=0 ; i< NpointsInFit ; i++) {
384  if( mvdhit[i])continue;
385  ii+= NStructRowsMax;
386 // fprintf(FMCS," SIGm%d OBJECT %g Am%d -1.\n SIGm%d Bm%d -1. Cm%d 1.\n SIGm%d Dm%d -1.\n",//---stampaggi
387 // i,1./Delta[i],i,i,i,i,i,i);//---stampaggi
388  Coefficients[ii]= 1./Delta[i];
389  Coefficients[ii+1]= -1.;
390  Coefficients[ii+2]= -1.;
391  Coefficients[ii+3]= 1.;
392  Coefficients[ii+4]= -1.;
393  }
394 
395 // Column variable DUM
396 //-----------stampaggi
397 /*
398  for(i=0 ; i< NpointsInFit ; i++) {
399  if( mvdhit[i]){
400  fprintf(FMCS," DUM A%d 1.\n DUM B%d 1.\n DUM C%d 1.\n",i,i,i);
401  fprintf(FMCS," DUM D%d 1.\n",i);
402  } else {
403  fprintf(FMCS," DUM Ap%d 1. Am%d 1.\n DUM Bp%d 1. Bm%d 1.\n DUM Cp%d 1. Cm%d 1\n",
404  i,i,i,i,i,i);
405  fprintf(FMCS," DUM Dp%d 1. Dm%d 1.\n",i,i);
406  }
407  }
408 */
409 //---fine stampaggi
410  ii+= NStructRowsMax;
411  for(i=0, iii=0 ; i< NpointsInFit ; i++) {
412  if( mvdhit[i]){
413  Coefficients[ii+iii]= 1.;
414  Coefficients[ii+iii+1]= 1.;
415  Coefficients[ii+iii+2]= 1.;
416  Coefficients[ii+iii+3]= 1.;
417  iii +=4;
418  } else {
419  Coefficients[ii+iii]= 1.;
420  Coefficients[ii+iii+1]= 1.;
421  Coefficients[ii+iii+2]= 1.;
422  Coefficients[ii+iii+3]= 1.;
423  Coefficients[ii+iii+4]= 1.;
424  Coefficients[ii+iii+5]= 1.;
425  Coefficients[ii+iii+6]= 1.;
426  Coefficients[ii+iii+7]= 1.;
427  iii +=8;
428  }
429  }
430 
431 //-----------
432 
433 //--- give the names to the Structural Variables
434 
435  sprintf(&auxStructVarName[0][0],"m1",i);
436  StructVarName[0] = &auxStructVarName[0][0];
437  sprintf(&auxStructVarName[1][0],"m2",i);
438  StructVarName[1] = &auxStructVarName[1][0];
439 
440  sprintf(&auxStructVarName[2][0],"q1",i);
441  StructVarName[2] = &auxStructVarName[2][0];
442 
443  sprintf(&auxStructVarName[3][0],"q2",i);
444  StructVarName[3] = &auxStructVarName[3][0];
445  for(i=0, ii=0; i< NpointsInFit ; i++) {
446  if( mvdhit[i]){
447  sprintf(&auxStructVarName[4+i][0],"lam%d",i);
448  StructVarName[4+i] = &auxStructVarName[4+i][0];
449 
450  sprintf(&auxStructVarName[4+nMvdHits+2*nSttHits+i][0],"SIG%d",i);
451  StructVarName[4+nMvdHits+2*nSttHits+i] = &auxStructVarName[4+nMvdHits+2*nSttHits+i][0];
452  } else {
453  sprintf(&auxStructVarName[4+i][0],"lamp%d",i);
454  StructVarName[4+i] = &auxStructVarName[4+i][0];
455 
456  sprintf(&auxStructVarName[4+NpointsInFit+ii][0],"lamm%d",i);
457  StructVarName[4+NpointsInFit+ii] = &auxStructVarName[4+NpointsInFit+ii][0];
458 
459  sprintf(&auxStructVarName[4+nMvdHits+2*nSttHits+i][0],"SIGp%d",i);
460  StructVarName[4+nMvdHits+2*nSttHits+i] = &auxStructVarName[4+nMvdHits+2*nSttHits+i][0];
461 
462  sprintf(&auxStructVarName[4+NpointsInFit+nMvdHits+2*nSttHits+ii][0],"SIGm%d",i);
463  StructVarName[4+NpointsInFit+nMvdHits+2*nSttHits+ii] =
464  &auxStructVarName[4+NpointsInFit+nMvdHits+2*nSttHits+ii][0];
465  ii++;
466  }
467  }
468 
469 
470  sprintf(&auxStructVarName[NStructVar-1][0],"DUM",i);
471  StructVarName[NStructVar-1] = &auxStructVarName[NStructVar-1][0];
472 
473 
474 
475 //--- give the names of those Rows in which the Structural Variables are present
476 
477 // for m1, m2, q1, q2
478  for(i=0; i< 4; i++){
479  for(j=0, ii=0; j< NpointsInFit;j++){
480  if( mvdhit[j]){
481  sprintf(&aux[i*NStructRowsMax+ii][0],"A%d",j);
482  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii]=&aux[i*NStructRowsMax+ii][0];
483  sprintf(&aux[i*NStructRowsMax+ii+1][0],"B%d",j);
484  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii+1]=&aux[i*NStructRowsMax+ii+1][0];
485  ii += 2;
486  } else {
487  sprintf(&aux[i*NStructRowsMax+ii][0],"Ap%d",j);
488  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii]=&aux[i*NStructRowsMax+ii][0];
489  sprintf(&aux[i*NStructRowsMax+ii+1][0],"Am%d",j);
490  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii+1]=&aux[i*NStructRowsMax+ii+1][0];
491  sprintf(&aux[i*NStructRowsMax+ii+2][0],"Bp%d",j);
492  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii+2]=&aux[i*NStructRowsMax+ii+2][0];
493  sprintf(&aux[i*NStructRowsMax+ii+3][0],"Bm%d",j);
494  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii+3]=&aux[i*NStructRowsMax+ii+3][0];
495  ii += 4;
496  }
497  }
498  }
499 
500 // now for the lamp* variables
501  for(i=0; i< NpointsInFit;i++){
502  if( mvdhit[i]){
503  sprintf(&aux[(i+4)*NStructRowsMax+0][0],"A%d",i);
504  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+0]= &aux[(i+4)*NStructRowsMax+0][0];
505  sprintf(&aux[(i+4)*NStructRowsMax+1][0],"B%d",i);
506  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+1]= &aux[(i+4)*NStructRowsMax+1][0];
507  sprintf(&aux[(i+4)*NStructRowsMax+2][0],"C%d",i);
508  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+2]= &aux[(i+4)*NStructRowsMax+2][0];
509  sprintf(&aux[(i+4)*NStructRowsMax+3][0],"D%d",i);
510  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+3]= &aux[(i+4)*NStructRowsMax+3][0];
511  } else {
512  sprintf(&aux[(i+4)*NStructRowsMax+0][0],"Ap%d",i);
513  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+0]= &aux[(i+4)*NStructRowsMax+0][0];
514  sprintf(&aux[(i+4)*NStructRowsMax+1][0],"Bp%d",i);
515  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+1]= &aux[(i+4)*NStructRowsMax+1][0];
516  sprintf(&aux[(i+4)*NStructRowsMax+2][0],"Cp%d",i);
517  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+2]= &aux[(i+4)*NStructRowsMax+2][0];
518  sprintf(&aux[(i+4)*NStructRowsMax+3][0],"Dp%d",i);
519  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+3]= &aux[(i+4)*NStructRowsMax+3][0];
520  sprintf(&aux[(i+4)*NStructRowsMax+4][0],"LM%d",i);
521  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+4]= &aux[(i+4)*NStructRowsMax+4][0];
522  }
523  }
524 
525 // now for the lamm* variables
526  for(i=0, ii=0; i< NpointsInFit;i++){
527  if( mvdhit[i]) continue;
528  sprintf(&aux[(ii+4+NpointsInFit)*NStructRowsMax][0],"Am%d",i);
529  NameRowsInWhichStructVarArePresent[(ii+4+NpointsInFit)*NStructRowsMax]=
530  &aux[(ii+4+NpointsInFit)*NStructRowsMax][0];
531  sprintf(&aux[(ii+4+NpointsInFit)*NStructRowsMax+1][0],"Bm%d",i);
532  NameRowsInWhichStructVarArePresent[(ii+4+NpointsInFit)*NStructRowsMax+1]=
533  &aux[(ii+4+NpointsInFit)*NStructRowsMax+1][0];
534  sprintf(&aux[(ii+4+NpointsInFit)*NStructRowsMax+2][0],"Cm%d",i);
535  NameRowsInWhichStructVarArePresent[(ii+4+NpointsInFit)*NStructRowsMax+2]=
536  &aux[(ii+4+NpointsInFit)*NStructRowsMax+2][0];
537  sprintf(&aux[(ii+4+NpointsInFit)*NStructRowsMax+3][0],"Dm%d",i);
538  NameRowsInWhichStructVarArePresent[(ii+4+NpointsInFit)*NStructRowsMax+3]=
539  &aux[(ii+4+NpointsInFit)*NStructRowsMax+3][0];
540  sprintf(&aux[(ii+4+NpointsInFit)*NStructRowsMax+4][0],"LM%d",i);
541  NameRowsInWhichStructVarArePresent[(ii+4+NpointsInFit)*NStructRowsMax+4]=
542  &aux[(ii+4+NpointsInFit)*NStructRowsMax+4][0];
543  ii++;
544  }
545 
546 // now for the SIGp* variables
547  for(i=0; i< NpointsInFit;i++){
548  if( mvdhit[i]) {
549 // sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax][0],"OBJECT");
550  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax]=
551  OBJECTname;
552  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1][0],"A%d",i);
553  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1]=
554  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1][0];
555  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2][0],"B%d",i);
556  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2]=
557  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2][0];
558  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3][0],"C%d",i);
559  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3]=
560  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3][0];
561  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4][0],"D%d",i);
562  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4]=
563  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4][0];
564  } else {
565 // sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax][0],"OBJECT");
566  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax]=
567  OBJECTname;
568  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1][0],"Ap%d",i);
569  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1]=
570  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1][0];
571  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2][0],"Bp%d",i);
572  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2]=
573  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2][0];
574  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3][0],"Cp%d",i);
575  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3]=
576  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3][0];
577  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4][0],"Dp%d",i);
578  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4]=
579  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4][0];
580  }
581  }
582 
583 // now for the SIGm* variables
584  for(i=0, ii=0; i< NpointsInFit;i++){
585  if( mvdhit[i]) continue;
586 // sprintf(&aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax][0],"OBJECT");
587  NameRowsInWhichStructVarArePresent[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax]=
588  OBJECTname;
589 
590  sprintf(&aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+1][0],"Am%d",i);
591  NameRowsInWhichStructVarArePresent[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+1]=
592  &aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+1][0];
593 
594  sprintf(&aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+2][0],"Bm%d",i);
595  NameRowsInWhichStructVarArePresent[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+2]=
596  &aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+2][0];
597 
598  sprintf(&aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+3][0],"Cm%d",i);
599  NameRowsInWhichStructVarArePresent[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+3]=
600  &aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+3][0];
601 
602  sprintf(&aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+4][0],"Dm%d",i);
603  NameRowsInWhichStructVarArePresent[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+4]=
604  &aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+4][0];
605  ii++;
606  }
607 
608 // now for the DUM variable
609  for(i=0, ii=0; i< NpointsInFit;i++){
610  if( mvdhit[i]) {
611  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii][0],"A%d",i);
612  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii]=
613  &aux[(NStructVar-1)*NStructRowsMax+ii][0];
614 
615  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+1][0],"B%d",i);
616  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+1]=
617  &aux[(NStructVar-1)*NStructRowsMax+ii+1][0];
618 
619  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+2][0],"C%d",i);
620  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+2]=
621  &aux[(NStructVar-1)*NStructRowsMax+ii+2][0];
622 
623  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+3][0],"D%d",i);
624  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+3]=
625  &aux[(NStructVar-1)*NStructRowsMax+ii+3][0];
626  ii += 4;
627  } else {
628  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii][0],"Ap%d",i);
629  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii]=
630  &aux[(NStructVar-1)*NStructRowsMax+ii][0];
631 
632  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+1][0],"Am%d",i);
633  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+1]=
634  &aux[(NStructVar-1)*NStructRowsMax+ii+1][0];
635 
636  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+2][0],"Bp%d",i);
637  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+2]=
638  &aux[(NStructVar-1)*NStructRowsMax+ii+2][0];
639 
640  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+3][0],"Bm%d",i);
641  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+3]=
642  &aux[(NStructVar-1)*NStructRowsMax+ii+3][0];
643  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+4][0],"Cp%d",i);
644  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+4]=
645  &aux[(NStructVar-1)*NStructRowsMax+ii+4][0];
646 
647  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+5][0],"Cm%d",i);
648  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+5]=
649  &aux[(NStructVar-1)*NStructRowsMax+ii+5][0];
650 
651  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+6][0],"Dp%d",i);
652  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+6]=
653  &aux[(NStructVar-1)*NStructRowsMax+ii+6][0];
654 
655  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+7][0],"Dm%d",i);
656  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+7]=
657  &aux[(NStructVar-1)*NStructRowsMax+ii+7][0];
658  ii += 8;
659  }
660  }
661 
662 
663 
664 //----------------- write the RHS section
665 
666 // fprintf(FMCS,"RHS\n");//---stampaggi
667  for(i=0, ii=0 ; i< NpointsInFit ; i++) {
668  if( mvdhit[i]){
669 //---stampaggi
670 /*
671  fprintf(FMCS," BOUND A%d %g B%d %g\n BOUND C%d %g D%d %g\n",
672  i, Oyy[i]+2.*M,i,
673  -Oyy[i]+2.*M,i,
674  Delta[i]+2.*M,i,M-Delta[i]+2.*M);
675 */
676 //---fine stampaggi
677  ValueB[ii] = Oyy[i]+2.*M;
678  ValueB[ii+1]= -Oyy[i]+2.*M;
679  ValueB[ii+2]= Delta[i]+2.*M;
680  ValueB[ii+3]= M-Delta[i]+2.*M;
681  ii += 4;
682  } else {
683 //---stampaggi
684 /*
685  fprintf(FMCS," BOUND Ap%d %g Bp%d %g\n BOUND Cp%d %g Dp%d %g\n",
686  i, Oyy[i]+DriftRadiusconformal[ i ]+2.*M,i,
687  -Oyy[i]-DriftRadiusconformal[ i ]+2.*M,i,
688  Delta[i]+2.*M,i,M-Delta[i]+2.*M);
689 */
690 //---fine stampaggi
691  ValueB[ii] = Oyy[i]+DriftRadiusconformal[ i ]+2.*M;
692  ValueB[ii+1]= -Oyy[i]-DriftRadiusconformal[ i ]+2.*M;
693  ValueB[ii+2]= Delta[i]+2.*M;
694  ValueB[ii+3]= M-Delta[i]+2.*M;
695 
696 
697 //---stampaggi
698 /*
699  fprintf(FMCS," BOUND Am%d %g Bm%d %g\n BOUND Cm%d %g Dm%d %g\n",
700  i, Oyy[i]-DriftRadiusconformal[ i ]+2.*M,i,
701  -Oyy[i]+DriftRadiusconformal[ i ]+2.*M,i,
702  Delta[i]+2.*M,i,M-Delta[i]+2.*M);
703  fprintf(FMCS," BOUND LM%d 1.\n",i);
704 */
705 //---fine stampaggi
706  ValueB[ii+4]= Oyy[i]-DriftRadiusconformal[ i ]+2.*M;
707  ValueB[ii+5]= -Oyy[i]+DriftRadiusconformal[ i ]+2.*M;
708  ValueB[ii+6]= Delta[i]+2.*M;
709  ValueB[ii+7]= M-Delta[i]+2.*M;
710  ValueB[ii+8]= 1.;
711  ii +=9;
712  }
713 
714  }
715 
716 
717 //----------------- write the RANGES section
718 
719 // fprintf(FMCS,"RANGES\n");//---stampaggi
720  for(i=0 , ii=0; i< NpointsInFit ; i++) {
721  if( mvdhit[i]) continue;
722 // fprintf(FMCS," RANGE LM%d 1.\n",i);//---stampaggi
723 //---
724  ValueRanges[ii]=1.;
725  sprintf(&auxNameRanges[ii][0],"LM%d",i);
726  NameRanges[ii]=&auxNameRanges[ii][0];
727  ii++;
728  }
729 
730 //----------------- write the BOUNDS section
731 // fprintf(FMCS,"BOUNDS\n");//---stampaggi
732 
733 
734  for(i=0 ; i< NpointsInFit ; i++) {
735  if( mvdhit[i]){
736 // fprintf(FMCS," BV Bounds lam%d\n", i);//---stampaggi
737  sprintf(&auxTypeofBound[i][0],"BV"); TypeofBound[i]= &auxTypeofBound[i][0];
738  sprintf(&auxBoundStructVarName[i][0],"lam%d",i);
739  } else {
740 // fprintf(FMCS," BV Bounds lamp%d\n", i);//---stampaggi
741  sprintf(&auxTypeofBound[i][0],"BV"); TypeofBound[i]= &auxTypeofBound[i][0];
742  sprintf(&auxBoundStructVarName[i][0],"lamp%d",i);
743  }
744 
745  BoundStructVarName[i]=&auxBoundStructVarName[i][0];
746  BoundValue[i]=0.;
747 
748  }
749 
750  for(i=0, ii=0 ; i< NpointsInFit ; i++) {
751  if( mvdhit[i]) continue;
752 // fprintf(FMCS," BV Bounds lamm%d\n", i);//---stampaggi
753  sprintf(&auxTypeofBound[ii+NpointsInFit][0],"BV");
754  TypeofBound[ii+NpointsInFit]= &auxTypeofBound[ii+NpointsInFit][0];
755  sprintf(&auxBoundStructVarName[ii+NpointsInFit][0],"lamm%d",i);
756  BoundStructVarName[ii+NpointsInFit]=&auxBoundStructVarName[ii+NpointsInFit][0];
757  BoundValue[ii+NpointsInFit]=0.;
758  ii++;
759  }
760 
761 // fprintf(FMCS," FX Bounds DUM %g\n",2.*M);//--stampaggi
762  sprintf(&auxTypeofBound[NpointsInFit+nSttHits][0],"FX");
763  TypeofBound[NpointsInFit+nSttHits]= &auxTypeofBound[NpointsInFit+nSttHits][0];
764 
765  sprintf(&auxTypeofBound[NpointsInFit+nSttHits][0],"FX");
766  TypeofBound[NpointsInFit+nSttHits]= &auxTypeofBound[NpointsInFit+nSttHits][0];
767 
768  sprintf(&auxBoundStructVarName[NpointsInFit+nSttHits][0],"DUM");
769  BoundStructVarName[NpointsInFit+nSttHits]=
770  &auxBoundStructVarName[NpointsInFit+nSttHits][0];
771  BoundValue[NpointsInFit+nSttHits]=2.*M;
772 //-----
773 
774 //------------------------------------------------------------stampaggi
775 /*
776  fprintf(FMCS,"ENDATA\n");
777  fclose(FMCS);
778 */
779 
780 if(istampa>=4){
781 
782 cout<<"n. punti nel fit "<<NpointsInFit<<endl;
783 
784 
785 cout<<"nRows "<<nRows<<endl;
786 for(int ic =0;ic<nRows; ic++){
787  cout<<"n. Row "<<ic<<", nameRows "<<nameRows[ic]<<", typeRows "<<typeRows[ic]<<endl;
788 }
789 
790 cout<<"NStructRowsMax, NStructVar "<<NStructRowsMax<<", "<<NStructVar<<endl;
791 for(int ic =0;ic<NStructVar; ic++){
792  cout<<"NRowsInWhichStructVarArePresent "<<NRowsInWhichStructVarArePresent[ic]
793  <<", nome var. strut. n."<<ic<<" = "
794  <<StructVarName[ic]<<endl;
795 
796  for(int jc=0; jc<NRowsInWhichStructVarArePresent[ic];jc++){
797  cout<<"n. "<<jc<<" NameRowsInWhichStructVarArePresent "
798  <<NameRowsInWhichStructVarArePresent[ic*NStructRowsMax+jc]<<endl;
799  }
800 }
801 
802 
803 cout<<"n Coefficient "<<21*nMvdHits+24*nSttHits<<endl;
804 iii=0;
805 for(int ic =0;ic<NStructVar; ic++){
806  cout<<"Struct. Var."<< StructVarName[ic] <<" e' presente in "
807  << NRowsInWhichStructVarArePresent[ic]<<" Rows;"<<endl;
808  for(ii=0;ii<NRowsInWhichStructVarArePresent[ic];ii++){
809 
810  cout<<"\tin Row "<<NameRowsInWhichStructVarArePresent[ic*NStructRowsMax+ii]
811  <<", ha Coefficient "<<Coefficients[ic*NStructRowsMax+ii]<<
812  " (n. sequenziale = "<<iii<<")"<<endl;
813  iii++;
814  }
815 }
816 
817 cout<<"n valuesB "<<nRows-1<<endl;
818 for(int ic =0;ic<nRows-1; ic++){
819  cout<<"n. "<<ic<<", valuesB "<<ValueB[ic]<<endl;
820 }
821 cout<<"n ranges "<<nRanges<<endl;
822 for(int ic =0;ic<nRanges; ic++){
823  cout<<"n. "<<ic<<", RANGES "<<ValueRanges[ic]<<endl;
824 }
825 cout<<"n Bounds "<<nBounds<<endl;
826 for(int ic =0;ic<nBounds; ic++){
827  cout<<"n. "<<ic<<", Bounds "<<BoundValue[ic]<<endl;
828  cout<<"n. "<<ic<<", Bound Type "<<TypeofBound[ic]<<endl;
829  cout<<"n. "<<ic<<", Bound Name "<<BoundStructVarName[ic]<<endl;
830 }
831  } // end of if(istampa
832 
833 //-------fine stampaggi
834 
835 //----------------------- funzioni chiamate direttamente
836 /*
837 cout<<"cavolo, da sttmvdtracking : nRows = "<<nRows<<", NStructVar = "<<
838  NStructVar<<", NStructRowsMax = "<<NStructRowsMax<<
839  ", NRowsInWhichStructVarArePresent = "<<
840  NRowsInWhichStructVarArePresent<<", nRanges = "<<nRanges
841  <<", nBounds = "<<nBounds<<endl;
842 */
843  int status= glp_main(
844  nRows,
845  nameRows,
846  typeRows, // ROWS info
847  NStructVar,
848  NStructRowsMax,
849  NRowsInWhichStructVarArePresent, // COLUMNS info
850  StructVarName,
851  NameRowsInWhichStructVarArePresent, // COLUMNS info
852  Coefficients, // COLUMNS info
853  ValueB, // RHS info
854  nRanges,
855  ValueRanges,
856  NameRanges, // RANGES info
857  nBounds,
858  BoundValue,
859  BoundStructVarName,
860  TypeofBound // BOUNDS info
861 // ,final_values,TIMEOUT
862  ,final_values
863  );
864 
865 
866  if(status != 0) { *Type=false; return -5 ;}
867 
868 
869 //--------stampaggi
870 if(istampa>=3){
871 printf("from FitHelixCylinder printout dopo glp_main -------------------------------\n");
872 printf(" number of structural variables %d\n",NStructVar);
873 int ica;
874 for(ica=0;ica<NStructVar;ica++){
875  printf("name of structural variable %s and its final value %g\n",
876  StructVarName[ica], final_values[ica]);
877 }
878 printf("from FitHelixCylinder printout dopo glp_main -------------------------------\n");
879 } // end of if(istampa
880 //--------fine stampaggi
881 
882 
883 
884 //----------------------- fine funzioni chiamate direttamente
885 
886 
887 /*
888  if(istampa>=3 && IVOLTE <= 20){
889  sprintf(stringa,
890 "/home/boca/panda/glpk/glpk-4.39/examples/glpsol --min -o soluztrack%dEvent%dstep%d GeneralParallelHitsConformeTraccia%dEvent%d.mcs",
891  0,IVOLTE,1,0,IVOLTE);
892  } else {
893  sprintf(stringa,
894 "/home/boca/panda/glpk/glpk-4.39/examples/glpsol --min -o soluztrack%dEvent%dstep%d GeneralParallelHitsConformeTraccia%dEvent%d.mcs >& /dev/null",
895  0, IVOLTE,1,0,IVOLTE,1);
896  }
897 */
898 
899  m1_result = final_values[0];
900  m2_result = final_values[1];
901  q1_result = final_values[2];
902  q2_result = final_values[3];
903 
904 if(istampa>2) cout<<"Results : m1 = "<<m1_result<<", m2= "<<m2_result<<", q1 = "<<q1_result<<
905  ", q2 = "<<q2_result<<endl;
906 
907 //--------- case in which the fit failed
908  if( final_values[0]==0. && final_values[1]==0. && final_values[2]==0. &&final_values[3]==0. )
909  { *Type = false; return -10;}
910 
911 //------------------------ transformation of the result in terms of ALFA, BETA, GAMMA
912 
913 
914  *qu = q1_result - q2_result;
915  *emme = m1_result-m2_result ;
916 
917  *pGamma = 0.;
918  if( fabs( *qu ) > 1.e-10) { // trajectory is a circle in XY space
919  *pAlfa = *emme/(*qu);
920  *pBeta = -1./(*qu);
921  *Type=true;
922 // now take into account the rotation and correct; the only affected quantities are ALFA and BETA
923  alfetta = *pAlfa;
924  *pAlfa = *pAlfa*cose - *pBeta*sine;
925  *pBeta = alfetta*sine + *pBeta*cose;
926 
927 
928  } else if(fabs(*emme)> 1.e-10) {// trajectory is a straight line in XY space of equation y= m*x
929  // the rotation first
930  angle = atan(*emme) + rotationangle;
931  if( fabs(cos(angle)) > 1.e-10 ) {
932  *pAlfa = 999999.;
933  *pBeta = -(*pAlfa)/tan(angle);
934  *Type=false;
935 
936  } else { // in this case the equation in XY plane is y = 0.
937  *pAlfa = 0.;
938  *pBeta = 999999.;
939  *Type=false;
940  }
941  } else { // in this case also the equation in XY plane is y = 0.
942  *pAlfa = 0.;
943  *pBeta = 999999.;
944  *Type=false;
945  } // end of if( fabs( *qu ) > 1.e-10)
946 
947 //------------------
948 
949 
950 // now take into account the displacement and correct
951  *pGamma += (trajectory_vertex[0]*trajectory_vertex[0]+ trajectory_vertex[1]*trajectory_vertex[1]
952  -*pAlfa*trajectory_vertex[0]-*pBeta*trajectory_vertex[1]);
953  *pAlfa -= 2.*trajectory_vertex[0];
954  *pBeta -= 2.*trajectory_vertex[1];
955 
956 
957  if(fabs(cose-*emme*sine)> 1.e-10) {
958  *qu=*qu/(cose-*emme*sine);
959  *emme=(*emme*cose+sine)/(cose-*emme*sine);
960  return 1;
961  } else { // in this case the equation is 0 = x+*qu .
962  if(fabs(sine+*emme*cose) < 1.e-10) {
963  cout<<" From FitHelixCylinder, equation of XY circle : X**2 + Y**2 =0,"
964  <<" situation impossible in principle! Returning -1"
965  <<endl;
966  return -1;
967  }
968 
969  *emme=1.;
970  *qu = *qu/(sine+*emme*cose);
971  return 99; // in this case the equation is 0 = x+*qu .
972  }
973 
974 
975 
976 
977 }
978 
979 
980 //----------end of function PndTrkGlpkFits::FitHelixCylinder
981 
982 
983 
984 //----------begin of function PndTrkGlpkFits::FitSZspace
985 
986 Short_t PndTrkGlpkFits::FitSZspace(
987  Short_t nSkewHitsinTrack,
988  Double_t *S,
989  Double_t *Z, //
990  Double_t *DriftRadius,
991  Double_t *ErrorDriftRadius, //
992  Double_t FInot,
993  Short_t NMAX,
994  Double_t *emme,
995  int IVOLTE
996  )
997 {
998 
999 
1000  // definition of variables for the glpsol solver
1001  // ROWS (for read_rows function)
1002  //
1003  Short_t NpointsInFit = nSkewHitsinTrack-NMAX <0 ? nSkewHitsinTrack : NMAX;
1004 
1005 
1006  bool mvdhit[NpointsInFit];
1007 
1008  Double_t ave,
1009  avex,
1010  avey,
1011  cose,
1012  sine,
1013  M = 50.,
1014  m_result,
1015  q_result,
1016  A,
1017  alfetta,
1018  angle,
1019  offsety,
1020  rotationangle,
1021  Oxx[NpointsInFit],
1022  Oyy[NpointsInFit],
1023  Delta[NpointsInFit];
1024 
1025  Short_t i, j, ii, iii, n, nSttHits, nMvdHits;
1026  Short_t Status;
1027 
1028  float m1_result,m2_result, q1_result,q2_result, A1_result, A2_result;
1029 
1030 // --
1031 
1032  if(nSkewHitsinTrack==0) {
1033  cout<<"from FitSZspace, Evento "<<IVOLTE<<endl;
1034  cout<<"from PndTrkGlpkFits::FitSZspace : no points in fit, return!\n";
1035  return -10;
1036  }
1037 
1038 
1039 /*
1040  ave=0.;
1041  avex=0.;
1042  avey=0.;
1043  n=0;
1044  for(i=0;i<nSkewHitsinTrack;i++){
1045  if( fabs(Z[ i ]) > 1.e-10&& DriftRadius[ i ]>0.){
1046  n++;
1047  ave += (S[ i ] - FInot)/Z[ i ];
1048  avex += Z[ i ];
1049  avey += (S[ i ] - FInot);
1050  }
1051  }
1052 
1053  if( n>0) {
1054  ave /=n;
1055  avex /=n;
1056  avey /=n;
1057  rotationangle = atan2(avey,avex);
1058  } else {
1059  rotationangle=PI/2.;
1060  }
1061 
1062 */
1063 
1064 
1065 
1066  rotationangle = PI/2.;
1067 
1068 
1069 // cose = cos(rotationangle);
1070 // sine = sin(rotationangle);
1071  cose = 0.;
1072  sine = 1.;
1073 
1074 
1075  // Delta[i] is actually used in the fit, and also Drift Radius.
1076 
1077  nSttHits = nMvdHits = 0;
1078  for(i=0;i<NpointsInFit; i++){
1079  Oxx[i] = Z[i]*cose +(S[i] - FInot)*sine;
1080  Oyy[i] = -Z[i]*sine +(S[i] - FInot)*cose;
1081  Delta[i] = ErrorDriftRadius[i];
1082 // Delta[i] = 2.*DriftRadius[i];
1083 
1084 
1085 
1086  if( DriftRadius[i]<0. )
1087  {
1088  mvdhit[i]=true;
1089  nMvdHits++;
1090  } else {
1091  mvdhit[i]=false;
1092  nSttHits++;
1093  }
1094  }
1095 
1096 
1097 
1098 
1099 
1100 //--------- calculation of # structural variables (see Gianluigi's logbook pag. 236) etc.etc.
1101 
1102  int NStructVar = 4 + 1 + nMvdHits * 2 + nSttHits *4 ;
1103 // m1,m2,q1,q2 DUM lam & SIG lamp & lamm & SIGp & SIGm
1104 
1105 
1106  int nRows = 1 + nMvdHits * 4 + nSttHits * 9;
1107 // OBJECT A,B,C,D Ap,Bp,Cp,Dp,Am,Bm,Cm,Dm,LM
1108 
1109 //---- creating the various service arrays
1110 
1111  char OBJECTname[8];
1112  sprintf(OBJECTname,"OBJECT");
1113 
1114  int typeRows[nRows];
1115  char * nameRows[nRows];
1116  char auxnameRows[nRows][5];
1117 
1118  int NStructRowsMax = 8*NpointsInFit ; // maximum number of ROWS in which a
1119  // structural variable (for instance M ) can be found
1120  double final_values[NStructVar];
1121  int NRowsInWhichStructVarArePresent[NStructVar];
1122  char *StructVarName[NStructVar];
1123  char auxStructVarName[NStructVar][8];
1124  char *NameRowsInWhichStructVarArePresent[NStructVar*NStructRowsMax];
1125  char aux[NStructVar*NStructRowsMax][5];
1126  double Coefficients[NStructVar*NStructRowsMax];
1127 
1128  //--------for RHS information
1129  double ValueB[nRows-1]; // -1 because OBJECT dowsn't have a RHS boundary.
1130  //--------for RANGES information
1131  int nRanges = nSttHits;
1132  double ValueRanges[nRanges];
1133  char *NameRanges[nRanges];
1134  char auxNameRanges[nRanges][20];
1135 
1136  //--------for BOUNDS information
1137  int nBounds=NpointsInFit+nSttHits+1+2;// +2 because q1 = q2 = 0 fixed.
1138  double BoundValue[nBounds];
1139  char *BoundStructVarName[nBounds];
1140  char auxBoundStructVarName[nBounds][8];
1141  char *TypeofBound[nBounds];
1142  char auxTypeofBound[nBounds][20];
1143  //--------end BOUNDS information
1144 
1145 //---------------------------------------------------
1146 
1147 
1148 
1149 
1150 //--- calculate array NRowsInWhichStructVarArePresent
1151  NRowsInWhichStructVarArePresent[0] = // this is for m1
1152  NRowsInWhichStructVarArePresent[1] = // this is for m2
1153  NRowsInWhichStructVarArePresent[2] = // this is for q1
1154  NRowsInWhichStructVarArePresent[3] = nMvdHits*2 + nSttHits *4; // this is for q2
1155  //--- the following is for the lam* (Mvd Hits) and lamp* (Stt hits) structural variables
1156  for(i=0,ii=0; i< NpointsInFit ; i++) {
1157  if( mvdhit[i]){
1158  NRowsInWhichStructVarArePresent[4+ii]= 4;
1159  } else {
1160  NRowsInWhichStructVarArePresent[4+ii]= 5;
1161  }
1162  ii++;
1163  }
1164 
1165  // the lamm* (Stt hits) if any.
1166  for(i=0; i< NpointsInFit ; i++) {
1167  if( !mvdhit[i]){
1168  NRowsInWhichStructVarArePresent[4+ii]= 5;
1169  ii++;
1170  }
1171  }
1172 
1173 
1174  //--- the following is for the SIG structural variables
1175  for(i=0 ; i< nMvdHits+2*nSttHits ; i++) {
1176  NRowsInWhichStructVarArePresent[4+ii+i]= 5;
1177  }
1178 //--- the following is for the DUM structural variable
1179  NRowsInWhichStructVarArePresent[4+ii+nMvdHits+2*nSttHits]= nMvdHits*4 + nSttHits*8;
1180 
1181 
1182 
1183 
1184 //----------------- write the ROWS section
1185 
1186 
1187 //-------stampaggi
1188 /*
1189  sprintf(nome,"GeneralSkewEvent%d.mcs", IVOLTE);
1190  FILE * FMCS = fopen(nome,"w");
1191  fprintf(FMCS,"NAME FIT\n");
1192 
1193 
1194  fprintf(FMCS,"ROWS\n");
1195  fprintf(FMCS," N OBJECT\n");
1196  for(i=0 ; i< NpointsInFit ; i++) {
1197  if( mvdhit[i]){
1198  fprintf(FMCS," L A%d\n L B%d\n L C%d\n L D%d\n",i,i,i,i);
1199  }else{
1200  fprintf(FMCS," L Ap%d\n L Bp%d\n L Cp%d\n L Dp%d\n",i,i,i,i);
1201  fprintf(FMCS," L Am%d\n L Bm%d\n L Cm%d\n L Dm%d\n G LM%d\n",i,i,i,i,i);
1202  }
1203  }
1204 */
1205 //--------------------------fine stampaggi
1206 
1207 
1208 
1209 
1210 // sprintf(&(auxnameRows[0][0]),"OBJECT");
1211 // nameRows[0]=&auxnameRows[0][0];
1212  nameRows[0]=OBJECTname;
1213  typeRows[0]=GLP_FR;
1214  for(i=0 , ii=0 ; i< NpointsInFit ; i++) {
1215 
1216  if( mvdhit[i]){
1217  typeRows[1+ii]=GLP_UP;typeRows[2+ii]=GLP_UP;typeRows[3+ii]=GLP_UP;typeRows[4+ii]=GLP_UP;
1218  typeRows[5+ii]=GLP_LO;
1219 
1220  sprintf(&(auxnameRows[1+ii][0]),"A%d",i); nameRows[1+ii]=&auxnameRows[1+ii][0];
1221  sprintf(&(auxnameRows[2+ii][0]),"B%d",i); nameRows[2+ii]=&auxnameRows[2+ii][0];
1222  sprintf(&(auxnameRows[3+ii][0]),"C%d",i); nameRows[3+ii]=&auxnameRows[3+ii][0];
1223  sprintf(&(auxnameRows[4+ii][0]),"D%d",i); nameRows[4+ii]=&auxnameRows[4+ii][0];
1224  ii+=4;
1225 
1226  } else {
1227  typeRows[1+ii]=GLP_UP;typeRows[2+ii]=GLP_UP;typeRows[3+ii]=GLP_UP;typeRows[4+ii]=GLP_UP;
1228  typeRows[5+ii]=GLP_UP;typeRows[6+ii]=GLP_UP;typeRows[7+ii]=GLP_UP;typeRows[8+ii]=GLP_UP;
1229  typeRows[9+ii]=GLP_LO;
1230 
1231  sprintf(&(auxnameRows[1+ii][0]),"Ap%d",i); nameRows[1+ii]=&auxnameRows[1+ii][0];
1232  sprintf(&(auxnameRows[2+ii][0]),"Bp%d",i); nameRows[2+ii]=&auxnameRows[2+ii][0];
1233  sprintf(&(auxnameRows[3+ii][0]),"Cp%d",i); nameRows[3+ii]=&auxnameRows[3+ii][0];
1234  sprintf(&(auxnameRows[4+ii][0]),"Dp%d",i); nameRows[4+ii]=&auxnameRows[4+ii][0];
1235  sprintf(&(auxnameRows[5+ii][0]),"Am%d",i); nameRows[5+ii]=&auxnameRows[5+ii][0];
1236  sprintf(&(auxnameRows[6+ii][0]),"Bm%d",i); nameRows[6+ii]=&auxnameRows[6+ii][0];
1237  sprintf(&(auxnameRows[7+ii][0]),"Cm%d",i); nameRows[7+ii]=&auxnameRows[7+ii][0];
1238  sprintf(&(auxnameRows[8+ii][0]),"Dm%d",i); nameRows[8+ii]=&auxnameRows[8+ii][0];
1239  sprintf(&(auxnameRows[9+ii][0]),"LM%d",i); nameRows[9+ii]=&auxnameRows[9+ii][0];
1240  ii+=9;
1241  }
1242  }
1243 
1244 
1245 
1246 
1247 //----------------- write the COLUMNS section
1248 
1249 // fprintf(FMCS,"COLUMNS\n"); //--------stampaggi
1250 
1251 // Column variable m1
1252 
1253  ii=0;
1254  for(i=0 ; i< NpointsInFit ; i++) {
1255 
1256  if( mvdhit[i]){
1257 //---stampaggi
1258 // fprintf(FMCS," m1 A%d %g\n m1 B%d %g\n",i,Oxx[i],i,-Oxx[i]);
1259 //-----stampaggi, fine
1260 
1261  Coefficients[ii]= Oxx[i];
1262  Coefficients[ii+1]= -Oxx[i];
1263  ii +=2;
1264  } else {
1265 //---stampaggi
1266 // fprintf(FMCS," m1 Ap%d %g Am%d %g\n m1 Bp%d %g Bm%d %g\n",
1267 // i,Oxx[i],i,Oxx[i],i,-Oxx[i],i,-Oxx[i]);
1268 //-----stampaggi, fine
1269 
1270  Coefficients[ii]= Oxx[i];
1271  Coefficients[ii+1]= Oxx[i];
1272  Coefficients[ii+2]= -Oxx[i];
1273  Coefficients[ii+3]= -Oxx[i];
1274  ii += 4;
1275  }
1276  }
1277 
1278 
1279 
1280 // Column variable m2
1281  for(i=0, ii=0; i< NpointsInFit ; i++) {
1282  if( mvdhit[i]){
1283 // fprintf(FMCS," m2 A%d %g\n m2 B%d %g\n",//---stampaggi
1284 // i,-Oxx[i],i,Oxx[i]);//---stampaggi
1285  Coefficients[NStructRowsMax+ii]= -Oxx[i];
1286  Coefficients[NStructRowsMax+ii+1]= Oxx[i];
1287  ii += 2;
1288  } else {
1289 // fprintf(FMCS," m2 Ap%d %g Am%d %g\n m2 Bp%d %g Bm%d %g\n",//---stampaggi
1290 // i,-Oxx[i],i,-Oxx[i],i,Oxx[i],i,Oxx[i]);//---stampaggi
1291  Coefficients[NStructRowsMax+ii]= -Oxx[i];
1292  Coefficients[NStructRowsMax+ii+1]= -Oxx[i];
1293  Coefficients[NStructRowsMax+ii+2]= Oxx[i];
1294  Coefficients[NStructRowsMax+ii+3]= Oxx[i];
1295  ii += 4;
1296  }
1297  }
1298 
1299 // Column variable q1
1300  for(i=0, ii=0 ; i< NpointsInFit ; i++) {
1301  if( mvdhit[i]){
1302 // fprintf(FMCS," q1 A%d 1.\n q1 B%d -1.\n",//---stampaggi
1303 // i,i);//---stampaggi
1304  Coefficients[2*NStructRowsMax+ii]= 1.;
1305  Coefficients[2*NStructRowsMax+ii+1]= -1.;
1306  ii +=2;
1307  } else {
1308 // fprintf(FMCS," q1 Ap%d 1. Am%d 1.\n q1 Bp%d -1. Bm%d -1.\n",//---stampaggi
1309 // i,i,i,i);//---stampaggi
1310  Coefficients[2*NStructRowsMax+ii]= 1.;
1311  Coefficients[2*NStructRowsMax+ii+1]= 1.;
1312  Coefficients[2*NStructRowsMax+ii+2]= -1.;
1313  Coefficients[2*NStructRowsMax+ii+3]= -1.;
1314  ii += 4;
1315  }
1316  }
1317 
1318 // Column variable q2
1319  for(i=0, ii=0 ; i< NpointsInFit ; i++) {
1320  if( mvdhit[i]){
1321 // fprintf(FMCS," q2 A%d -1.\n q2 B%d 1.\n",//---stampaggi
1322 // i,i);//---stampaggi
1323  Coefficients[3*NStructRowsMax+ii]= -1.;
1324  Coefficients[3*NStructRowsMax+ii+1]= 1.;
1325  ii += 2;
1326  } else {
1327 // fprintf(FMCS," q2 Ap%d -1. Am%d -1.\n q2 Bp%d 1. Bm%d 1.\n",//---stampaggi
1328 // i,i,i,i);//---stampaggi
1329  Coefficients[3*NStructRowsMax+ii]= -1.;
1330  Coefficients[3*NStructRowsMax+ii+1]= -1.;
1331  Coefficients[3*NStructRowsMax+ii+2]= 1.;
1332  Coefficients[3*NStructRowsMax+ii+3]= 1.;
1333  ii += 4;
1334  }
1335 
1336  }
1337 
1338 // Column variable lambdap(i)
1339  for(i=0 ; i< NpointsInFit ; i++) {
1340  ii=(4+i)*NStructRowsMax;
1341  Coefficients[ii]= -M;
1342  Coefficients[ii+1]= -M;
1343  Coefficients[ii+2]= -M;
1344  Coefficients[ii+3]= M;
1345 /*
1346  if( mvdhit[i]){
1347  fprintf(FMCS," lam%d A%d %g B%d %g\n lam%d C%d %g D%d %g\n",//---stampaggi
1348  i,i,-M,i,-M, i , i,-M, i, M);//---stampaggi
1349  } else {
1350  fprintf(FMCS," lamp%d Ap%d %g Bp%d %g\n lamp%d Cp%d %g Dp%d %g\n lamp%d LM%d 1.\n",//---stampaggi
1351  i,i,-M,i,-M, i , i,-M, i, M, i,i);//---stampaggi
1352  }
1353 */
1354  if(! mvdhit[i]) Coefficients[ii+4]= 1.;
1355 
1356  }
1357 
1358 
1359 
1360 // Column variable lambdam(i)
1361  for(i=0 ; i< NpointsInFit ; i++) {
1362  if( mvdhit[i]) continue;
1363  ii+= NStructRowsMax;
1364 // fprintf(FMCS," lamm%d Am%d %g Bm%d %g\n lamm%d Cm%d %g Dm%d %g\n lamm%d LM%d 1.\n",//---stampaggi
1365 // i,i,-M,i,-M, i,i , -M, i, M, i,i);//---stampaggi
1366  Coefficients[ii]= -M;
1367  Coefficients[ii+1]= -M;
1368  Coefficients[ii+2]= -M;
1369  Coefficients[ii+3]= M;
1370  Coefficients[ii+4]= 1.;
1371  }
1372 // Column variable SIGp(i)
1373  for(i=0; i< NpointsInFit ; i++) {
1374  ii+= NStructRowsMax;
1375 
1376 /*
1377  if( mvdhit[i]){
1378  fprintf(FMCS," SIG%d OBJECT %g A%d -1.\n SIG%d B%d -1. C%d 1.\n SIG%d D%d -1.\n",//---stampaggi
1379  i,1./Delta[i],i,i,i,i,i,i);//---stampaggi
1380 
1381  } else {
1382  fprintf(FMCS," SIGp%d OBJECT %g Ap%d -1.\n SIGp%d Bp%d -1. Cp%d 1.\n SIGp%d Dp%d -1.\n",//---stampaggi
1383  i,1./Delta[i],i,i,i,i,i,i);//---stampaggi
1384  }
1385 */
1386 
1387 
1388  Coefficients[ii]= 1./Delta[i];
1389  Coefficients[ii+1]= -1.;
1390  Coefficients[ii+2]= -1.;
1391  Coefficients[ii+3]= 1.;
1392  Coefficients[ii+4]= -1.;
1393 
1394 
1395  }
1396 // Column variable SIGm(i)
1397  for(i=0 ; i< NpointsInFit ; i++) {
1398  if( mvdhit[i])continue;
1399  ii+= NStructRowsMax;
1400 /*
1401  fprintf(FMCS," SIGm%d OBJECT %g Am%d -1.\n SIGm%d Bm%d -1. Cm%d 1.\n SIGm%d Dm%d -1.\n",//---stampaggi
1402  i,1./Delta[i],i,i,i,i,i,i);//---stampaggi
1403 */
1404  Coefficients[ii]= 1./Delta[i];
1405  Coefficients[ii+1]= -1.;
1406  Coefficients[ii+2]= -1.;
1407  Coefficients[ii+3]= 1.;
1408  Coefficients[ii+4]= -1.;
1409  }
1410 
1411 // Column variable DUM
1412 //-----------stampaggi
1413 /*
1414  for(i=0 ; i< NpointsInFit ; i++) {
1415  if( mvdhit[i]){
1416  fprintf(FMCS," DUM A%d 1.\n DUM B%d 1.\n DUM C%d 1.\n",i,i,i);
1417  fprintf(FMCS," DUM D%d 1.\n",i);
1418  } else {
1419  fprintf(FMCS," DUM Ap%d 1. Am%d 1.\n DUM Bp%d 1. Bm%d 1.\n DUM Cp%d 1. Cm%d 1\n",
1420  i,i,i,i,i,i);
1421  fprintf(FMCS," DUM Dp%d 1. Dm%d 1.\n",i,i);
1422  }
1423  }
1424 */
1425 //---fine stampaggi
1426  ii+= NStructRowsMax;
1427  for(i=0, iii=0 ; i< NpointsInFit ; i++) {
1428  if( mvdhit[i]){
1429  Coefficients[ii+iii]= 1.;
1430  Coefficients[ii+iii+1]= 1.;
1431  Coefficients[ii+iii+2]= 1.;
1432  Coefficients[ii+iii+3]= 1.;
1433  iii +=4;
1434  } else {
1435  Coefficients[ii+iii]= 1.;
1436  Coefficients[ii+iii+1]= 1.;
1437  Coefficients[ii+iii+2]= 1.;
1438  Coefficients[ii+iii+3]= 1.;
1439  Coefficients[ii+iii+4]= 1.;
1440  Coefficients[ii+iii+5]= 1.;
1441  Coefficients[ii+iii+6]= 1.;
1442  Coefficients[ii+iii+7]= 1.;
1443  iii +=8;
1444  }
1445  }
1446 
1447 //-----------
1448 
1449 //--- give the names to the Structural Variables
1450 
1451  sprintf(&auxStructVarName[0][0],"m1",i);
1452  StructVarName[0] = &auxStructVarName[0][0];
1453  sprintf(&auxStructVarName[1][0],"m2",i);
1454  StructVarName[1] = &auxStructVarName[1][0];
1455 
1456  sprintf(&auxStructVarName[2][0],"q1",i);
1457  StructVarName[2] = &auxStructVarName[2][0];
1458 
1459  sprintf(&auxStructVarName[3][0],"q2",i);
1460  StructVarName[3] = &auxStructVarName[3][0];
1461  for(i=0, ii=0; i< NpointsInFit ; i++) {
1462  if( mvdhit[i]){
1463  sprintf(&auxStructVarName[4+i][0],"lam%d",i);
1464  StructVarName[4+i] = &auxStructVarName[4+i][0];
1465 
1466  sprintf(&auxStructVarName[4+nMvdHits+2*nSttHits+i][0],"SIG%d",i);
1467  StructVarName[4+nMvdHits+2*nSttHits+i] = &auxStructVarName[4+nMvdHits+2*nSttHits+i][0];
1468  } else {
1469  sprintf(&auxStructVarName[4+i][0],"lamp%d",i);
1470  StructVarName[4+i] = &auxStructVarName[4+i][0];
1471 
1472  sprintf(&auxStructVarName[4+NpointsInFit+ii][0],"lamm%d",i);
1473  StructVarName[4+NpointsInFit+ii] = &auxStructVarName[4+NpointsInFit+ii][0];
1474 
1475  sprintf(&auxStructVarName[4+nMvdHits+2*nSttHits+i][0],"SIGp%d",i);
1476  StructVarName[4+nMvdHits+2*nSttHits+i] = &auxStructVarName[4+nMvdHits+2*nSttHits+i][0];
1477 
1478  sprintf(&auxStructVarName[4+NpointsInFit+nMvdHits+2*nSttHits+ii][0],"SIGm%d",i);
1479  StructVarName[4+NpointsInFit+nMvdHits+2*nSttHits+ii] =
1480  &auxStructVarName[4+NpointsInFit+nMvdHits+2*nSttHits+ii][0];
1481  ii++;
1482  }
1483  }
1484 
1485 
1486  sprintf(&auxStructVarName[NStructVar-1][0],"DUM",i);
1487  StructVarName[NStructVar-1] = &auxStructVarName[NStructVar-1][0];
1488 
1489 
1490 
1491 //--- give the names of those Rows in which the Structural Variables are present
1492 
1493 // for m1, m2, q1, q2
1494  for(i=0; i< 4; i++){
1495  for(j=0, ii=0; j< NpointsInFit;j++){
1496  if( mvdhit[j]){
1497  sprintf(&aux[i*NStructRowsMax+ii][0],"A%d",j);
1498  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii]=&aux[i*NStructRowsMax+ii][0];
1499  sprintf(&aux[i*NStructRowsMax+ii+1][0],"B%d",j);
1500  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii+1]=&aux[i*NStructRowsMax+ii+1][0];
1501  ii += 2;
1502  } else {
1503  sprintf(&aux[i*NStructRowsMax+ii][0],"Ap%d",j);
1504  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii]=&aux[i*NStructRowsMax+ii][0];
1505  sprintf(&aux[i*NStructRowsMax+ii+1][0],"Am%d",j);
1506  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii+1]=&aux[i*NStructRowsMax+ii+1][0];
1507  sprintf(&aux[i*NStructRowsMax+ii+2][0],"Bp%d",j);
1508  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii+2]=&aux[i*NStructRowsMax+ii+2][0];
1509  sprintf(&aux[i*NStructRowsMax+ii+3][0],"Bm%d",j);
1510  NameRowsInWhichStructVarArePresent[i*NStructRowsMax+ii+3]=&aux[i*NStructRowsMax+ii+3][0];
1511  ii += 4;
1512  }
1513  }
1514  }
1515 
1516 // now for the lamp* variables
1517  for(i=0; i< NpointsInFit;i++){
1518  if( mvdhit[i]){
1519  sprintf(&aux[(i+4)*NStructRowsMax+0][0],"A%d",i);
1520  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+0]= &aux[(i+4)*NStructRowsMax+0][0];
1521  sprintf(&aux[(i+4)*NStructRowsMax+1][0],"B%d",i);
1522  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+1]= &aux[(i+4)*NStructRowsMax+1][0];
1523  sprintf(&aux[(i+4)*NStructRowsMax+2][0],"C%d",i);
1524  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+2]= &aux[(i+4)*NStructRowsMax+2][0];
1525  sprintf(&aux[(i+4)*NStructRowsMax+3][0],"D%d",i);
1526  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+3]= &aux[(i+4)*NStructRowsMax+3][0];
1527  } else {
1528  sprintf(&aux[(i+4)*NStructRowsMax+0][0],"Ap%d",i);
1529  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+0]= &aux[(i+4)*NStructRowsMax+0][0];
1530  sprintf(&aux[(i+4)*NStructRowsMax+1][0],"Bp%d",i);
1531  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+1]= &aux[(i+4)*NStructRowsMax+1][0];
1532  sprintf(&aux[(i+4)*NStructRowsMax+2][0],"Cp%d",i);
1533  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+2]= &aux[(i+4)*NStructRowsMax+2][0];
1534  sprintf(&aux[(i+4)*NStructRowsMax+3][0],"Dp%d",i);
1535  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+3]= &aux[(i+4)*NStructRowsMax+3][0];
1536  sprintf(&aux[(i+4)*NStructRowsMax+4][0],"LM%d",i);
1537  NameRowsInWhichStructVarArePresent[(i+4)*NStructRowsMax+4]= &aux[(i+4)*NStructRowsMax+4][0];
1538  }
1539  }
1540 
1541 // now for the lamm* variables
1542  for(i=0, ii=0; i< NpointsInFit;i++){
1543  if( mvdhit[i]) continue;
1544  sprintf(&aux[(ii+4+NpointsInFit)*NStructRowsMax][0],"Am%d",i);
1545  NameRowsInWhichStructVarArePresent[(ii+4+NpointsInFit)*NStructRowsMax]= &aux[(ii+4+NpointsInFit)*NStructRowsMax][0];
1546  sprintf(&aux[(ii+4+NpointsInFit)*NStructRowsMax+1][0],"Bm%d",i);
1547  NameRowsInWhichStructVarArePresent[(ii+4+NpointsInFit)*NStructRowsMax+1]= &aux[(ii+4+NpointsInFit)*NStructRowsMax+1][0];
1548  sprintf(&aux[(ii+4+NpointsInFit)*NStructRowsMax+2][0],"Cm%d",i);
1549  NameRowsInWhichStructVarArePresent[(ii+4+NpointsInFit)*NStructRowsMax+2]= &aux[(ii+4+NpointsInFit)*NStructRowsMax+2][0];
1550  sprintf(&aux[(ii+4+NpointsInFit)*NStructRowsMax+3][0],"Dm%d",i);
1551  NameRowsInWhichStructVarArePresent[(ii+4+NpointsInFit)*NStructRowsMax+3]= &aux[(ii+4+NpointsInFit)*NStructRowsMax+3][0];
1552  sprintf(&aux[(ii+4+NpointsInFit)*NStructRowsMax+4][0],"LM%d",i);
1553  NameRowsInWhichStructVarArePresent[(ii+4+NpointsInFit)*NStructRowsMax+4]= &aux[(ii+4+NpointsInFit)*NStructRowsMax+4][0];
1554  ii++;
1555  }
1556 
1557 // now for the SIGp* variables
1558  for(i=0; i< NpointsInFit;i++){
1559  if( mvdhit[i]) {
1560 // sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax][0],"OBJECT");
1561  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax]=
1562  OBJECTname;
1563  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1][0],"A%d",i);
1564  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1]=
1565  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1][0];
1566  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2][0],"B%d",i);
1567  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2]=
1568  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2][0];
1569  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3][0],"C%d",i);
1570  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3]=
1571  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3][0];
1572  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4][0],"D%d",i);
1573  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4]=
1574  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4][0];
1575  } else {
1576 // sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax][0],"OBJECT");
1577  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax]=
1578  OBJECTname;
1579  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1][0],"Ap%d",i);
1580  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1]=
1581  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+1][0];
1582  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2][0],"Bp%d",i);
1583  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2]=
1584  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+2][0];
1585  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3][0],"Cp%d",i);
1586  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3]=
1587  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+3][0];
1588  sprintf(&aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4][0],"Dp%d",i);
1589  NameRowsInWhichStructVarArePresent[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4]=
1590  &aux[(i+4+NpointsInFit+nSttHits)*NStructRowsMax+4][0];
1591  }
1592  }
1593 
1594 // now for the SIGm* variables
1595  for(i=0, ii=0; i< NpointsInFit;i++){
1596  if( mvdhit[i]) continue;
1597 // sprintf(&aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax][0],"OBJECT");
1598  NameRowsInWhichStructVarArePresent[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax]=
1599  OBJECTname;
1600 
1601  sprintf(&aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+1][0],"Am%d",i);
1602  NameRowsInWhichStructVarArePresent[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+1]=
1603  &aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+1][0];
1604 
1605  sprintf(&aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+2][0],"Bm%d",i);
1606  NameRowsInWhichStructVarArePresent[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+2]=
1607  &aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+2][0];
1608 
1609  sprintf(&aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+3][0],"Cm%d",i);
1610  NameRowsInWhichStructVarArePresent[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+3]=
1611  &aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+3][0];
1612 
1613  sprintf(&aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+4][0],"Dm%d",i);
1614  NameRowsInWhichStructVarArePresent[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+4]=
1615  &aux[(ii+4+2*NpointsInFit+nSttHits)*NStructRowsMax+4][0];
1616  ii++;
1617  }
1618 
1619 // now for the DUM variable
1620  for(i=0, ii=0; i< NpointsInFit;i++){
1621  if( mvdhit[i]) {
1622  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii][0],"A%d",i);
1623  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii]=
1624  &aux[(NStructVar-1)*NStructRowsMax+ii][0];
1625 
1626  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+1][0],"B%d",i);
1627  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+1]=
1628  &aux[(NStructVar-1)*NStructRowsMax+ii+1][0];
1629 
1630  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+2][0],"C%d",i);
1631  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+2]=
1632  &aux[(NStructVar-1)*NStructRowsMax+ii+2][0];
1633 
1634  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+3][0],"D%d",i);
1635  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+3]=
1636  &aux[(NStructVar-1)*NStructRowsMax+ii+3][0];
1637  ii += 4;
1638  } else {
1639  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii][0],"Ap%d",i);
1640  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii]=
1641  &aux[(NStructVar-1)*NStructRowsMax+ii][0];
1642 
1643  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+1][0],"Am%d",i);
1644  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+1]=
1645  &aux[(NStructVar-1)*NStructRowsMax+ii+1][0];
1646 
1647  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+2][0],"Bp%d",i);
1648  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+2]=
1649  &aux[(NStructVar-1)*NStructRowsMax+ii+2][0];
1650 
1651  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+3][0],"Bm%d",i);
1652  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+3]=
1653  &aux[(NStructVar-1)*NStructRowsMax+ii+3][0];
1654  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+4][0],"Cp%d",i);
1655  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+4]=
1656  &aux[(NStructVar-1)*NStructRowsMax+ii+4][0];
1657 
1658  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+5][0],"Cm%d",i);
1659  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+5]=
1660  &aux[(NStructVar-1)*NStructRowsMax+ii+5][0];
1661 
1662  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+6][0],"Dp%d",i);
1663  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+6]=
1664  &aux[(NStructVar-1)*NStructRowsMax+ii+6][0];
1665 
1666  sprintf(&aux[(NStructVar-1)*NStructRowsMax+ii+7][0],"Dm%d",i);
1667  NameRowsInWhichStructVarArePresent[(NStructVar-1)*NStructRowsMax+ii+7]=
1668  &aux[(NStructVar-1)*NStructRowsMax+ii+7][0];
1669  ii += 8;
1670  }
1671  }
1672 
1673 
1674 
1675 //----------------- write the RHS section
1676 
1677 // fprintf(FMCS,"RHS\n");//---stampaggi
1678  for(i=0, ii=0 ; i< NpointsInFit ; i++) {
1679  if( mvdhit[i]){
1680 //---stampaggi
1681 /*
1682  fprintf(FMCS," BOUND A%d %g B%d %g\n BOUND C%d %g D%d %g\n",
1683  i, Oyy[i]+2.*M,i,
1684  -Oyy[i]+2.*M,i,
1685  Delta[i]+2.*M,i,M-Delta[i]+2.*M);
1686 */
1687 //---fine stampaggi
1688  ValueB[ii] = Oyy[i]+2.*M;
1689  ValueB[ii+1]= -Oyy[i]+2.*M;
1690  ValueB[ii+2]= Delta[i]+2.*M;
1691  ValueB[ii+3]= M-Delta[i]+2.*M;
1692  ii += 4;
1693  } else {
1694 //---stampaggi
1695 /*
1696  fprintf(FMCS," BOUND Ap%d %g Bp%d %g\n BOUND Cp%d %g Dp%d %g\n",
1697  i, Oyy[i]+DriftRadius[ i ]+2.*M,i,
1698  -Oyy[i]-DriftRadius[ i ]+2.*M,i,
1699  Delta[i]+2.*M,i,M-Delta[i]+2.*M);
1700 */
1701 //---fine stampaggi
1702  ValueB[ii] = Oyy[i]+DriftRadius[ i ]+2.*M;
1703  ValueB[ii+1]= -Oyy[i]-DriftRadius[ i ]+2.*M;
1704  ValueB[ii+2]= Delta[i]+2.*M;
1705  ValueB[ii+3]= M-Delta[i]+2.*M;
1706 
1707 
1708 //---stampaggi
1709 /*
1710  fprintf(FMCS," BOUND Am%d %g Bm%d %g\n BOUND Cm%d %g Dm%d %g\n",
1711  i, Oyy[i]-DriftRadius[ i ]+2.*M,i,
1712  -Oyy[i]+DriftRadius[ i ]+2.*M,i,
1713  Delta[i]+2.*M,i,M-Delta[i]+2.*M);f
1714  fprintf(FMCS," BOUND LM%d 1.\n",i);
1715 */
1716 //---fine stampaggi
1717  ValueB[ii+4]= Oyy[i]-DriftRadius[ i ]+2.*M;
1718  ValueB[ii+5]= -Oyy[i]+DriftRadius[ i ]+2.*M;
1719  ValueB[ii+6]= Delta[i]+2.*M;
1720  ValueB[ii+7]= M-Delta[i]+2.*M;
1721  ValueB[ii+8]= 1.;
1722  ii +=9;
1723  }
1724 
1725  }
1726 
1727 
1728 //----------------- write the RANGES section
1729 
1730 // fprintf(FMCS,"RANGES\n");//---stampaggi
1731  for(i=0 , ii=0; i< NpointsInFit ; i++) {
1732  if( mvdhit[i]) continue;
1733 // fprintf(FMCS," RANGE LM%d 1.\n",i);//---stampaggi
1734 //---
1735  ValueRanges[ii]=1.;
1736  sprintf(&auxNameRanges[ii][0],"LM%d",i);
1737  NameRanges[ii]=&auxNameRanges[ii][0];
1738  ii++;
1739  }
1740 
1741 //----------------- write the BOUNDS section
1742 // fprintf(FMCS,"BOUNDS\n");//---stampaggi
1743 
1744 
1745  for(i=0 ; i< NpointsInFit ; i++) {
1746  if( mvdhit[i]){
1747 // fprintf(FMCS," BV Bounds lam%d\n", i);//---stampaggi
1748  sprintf(&auxTypeofBound[i][0],"BV"); TypeofBound[i]= &auxTypeofBound[i][0];
1749  sprintf(&auxBoundStructVarName[i][0],"lam%d",i);
1750  } else {
1751 // fprintf(FMCS," BV Bounds lamp%d\n", i);//---stampaggi
1752  sprintf(&auxTypeofBound[i][0],"BV"); TypeofBound[i]= &auxTypeofBound[i][0];
1753  sprintf(&auxBoundStructVarName[i][0],"lamp%d",i);
1754  }
1755 
1756  BoundStructVarName[i]=&auxBoundStructVarName[i][0];
1757  BoundValue[i]=0.;
1758 
1759  }
1760 
1761  for(i=0, ii=0 ; i< NpointsInFit ; i++) {
1762  if( mvdhit[i]) continue;
1763 // fprintf(FMCS," BV Bounds lamm%d\n", i);//---stampaggi
1764  sprintf(&auxTypeofBound[ii+NpointsInFit][0],"BV");
1765  TypeofBound[ii+NpointsInFit]= &auxTypeofBound[ii+NpointsInFit][0];
1766  sprintf(&auxBoundStructVarName[ii+NpointsInFit][0],"lamm%d",i);
1767  BoundStructVarName[ii+NpointsInFit]=&auxBoundStructVarName[ii+NpointsInFit][0];
1768  BoundValue[ii+NpointsInFit]=0.;
1769  ii++;
1770  }
1771 
1772 // fprintf(FMCS," FX Bounds DUM %g\n",2.*M);//--stampaggi
1773  sprintf(&auxTypeofBound[NpointsInFit+nSttHits][0],"FX");
1774  TypeofBound[NpointsInFit+nSttHits]= &auxTypeofBound[NpointsInFit+nSttHits][0];
1775 
1776  sprintf(&auxTypeofBound[NpointsInFit+nSttHits][0],"FX");
1777  TypeofBound[NpointsInFit+nSttHits]= &auxTypeofBound[NpointsInFit+nSttHits][0];
1778 
1779  sprintf(&auxBoundStructVarName[NpointsInFit+nSttHits][0],"DUM");
1780  BoundStructVarName[NpointsInFit+nSttHits]=&auxBoundStructVarName[NpointsInFit+nSttHits][0];
1781  BoundValue[NpointsInFit+nSttHits]=2.*M;
1782 
1783 // fixing q1
1784 // fprintf(FMCS," FX Bounds q1 %g\n",0.);//--stampaggi
1785  sprintf(&auxTypeofBound[NpointsInFit+nSttHits+1][0],"FX");
1786  TypeofBound[NpointsInFit+nSttHits+1]=&auxTypeofBound[NpointsInFit+nSttHits+1][0];
1787  sprintf(&auxBoundStructVarName[NpointsInFit+nSttHits+1][0],"q1");
1788  BoundStructVarName[NpointsInFit+nSttHits+1]=&auxBoundStructVarName[NpointsInFit+nSttHits+1][0];
1789  BoundValue[NpointsInFit+nSttHits+1]= 0.;
1790 // fixing q2
1791 // fprintf(FMCS," FX Bounds q2 %g\n",0.);//--stampaggi
1792  sprintf(&auxTypeofBound[NpointsInFit+nSttHits+2][0],"FX");
1793  TypeofBound[NpointsInFit+nSttHits+2]=&auxTypeofBound[NpointsInFit+nSttHits+2][0];
1794  sprintf(&auxBoundStructVarName[NpointsInFit+nSttHits+2][0],"q2");
1795  BoundStructVarName[NpointsInFit+nSttHits+2]=&auxBoundStructVarName[NpointsInFit+nSttHits+2][0];
1796  BoundValue[NpointsInFit+nSttHits+2]= 0.;
1797 
1798 //-----
1799 
1800 //------------------------------------------------------------stampaggi
1801 /*
1802  fprintf(FMCS,"ENDATA\n");
1803  fclose(FMCS);
1804 */
1805 
1806 /*
1807 cout<<"n. punti nel fit "<<NpointsInFit<<endl;
1808 
1809 
1810 cout<<"nRows "<<nRows<<endl;
1811 for(int ic =0;ic<nRows; ic++){
1812  cout<<"n. Row "<<ic<<", nameRows "<<nameRows[ic]<<", typeRows "<<typeRows[ic]<<endl;
1813 }
1814 
1815 cout<<"NStructRowsMax = "<<NStructRowsMax<<endl;
1816 cout<<"NStructVar "<<NStructVar<<" e loro elenco "<<endl;
1817 for(int ic =0;ic<NStructVar; ic++){
1818  cout<<"\tvar. n. "<<ic<<", nome = "<<StructVarName[ic]<<endl;
1819 }
1820 
1821 
1822 
1823 for(int ic =0;ic<NStructVar; ic++){
1824  cout<<"NRowsInWhichStructVarArePresent "<<NRowsInWhichStructVarArePresent[ic]
1825  <<", nome var. strut. n."<<ic<<" = "
1826  <<StructVarName[ic]<<endl;
1827 
1828  for(int jc=0; jc<NRowsInWhichStructVarArePresent[ic];jc++){
1829  cout<<"n. "<<jc<<" NameRowsInWhichStructVarArePresent "<<NameRowsInWhichStructVarArePresent[ic*NStructRowsMax+jc]<<endl;
1830  }
1831 }
1832 
1833 
1834 cout<<"n Coefficient "<<21*nMvdHits+24*nSttHits<<endl;
1835 iii=0;
1836 for(int ic =0;ic<NStructVar; ic++){
1837  cout<<"Struct. Var."<< StructVarName[ic] <<" e' presente in "<< NRowsInWhichStructVarArePresent[ic]
1838  <<" Rows;"<<endl;
1839  for(ii=0;ii<NRowsInWhichStructVarArePresent[ic];ii++){
1840 
1841  cout<<"\tin Row "<<NameRowsInWhichStructVarArePresent[ic*NStructRowsMax+ii]
1842  <<", ha Coefficient "<<Coefficients[ic*NStructRowsMax+ii]<<
1843  " (n. sequenziale = "<<iii<<")"<<endl;
1844  iii++;
1845  }
1846 }
1847 
1848 cout<<"n valuesB "<<nRows-1<<endl;
1849 for(int ic =0;ic<nRows-1; ic++){
1850  cout<<"n. "<<ic<<", valuesB "<<ValueB[ic]<<endl;
1851 }
1852 cout<<"n ranges "<<nRanges<<endl;
1853 for(int ic =0;ic<nRanges; ic++){
1854  cout<<"n. "<<ic<<", RANGES "<<ValueRanges[ic]<<endl;
1855 }
1856 cout<<"n Bounds "<<nBounds<<endl;
1857 for(int ic =0;ic<nBounds; ic++){
1858  cout<<"n. "<<ic<<", Bounds "<<BoundValue[ic]<<endl;
1859  cout<<"n. "<<ic<<", Bound Type "<<TypeofBound[ic]<<endl;
1860  cout<<"n. "<<ic<<", Bound Name "<<BoundStructVarName[ic]<<endl;
1861 }
1862 */
1863 
1864 //-------fine stampaggi
1865 
1866 //----------------------- funzioni chiamate direttamente
1867 /*
1868 cout<<"cavolo2, da sttmvdtracking : nRows = "<<nRows<<", NStructVar = "<<
1869  NStructVar<<", NStructRowsMax = "<<NStructRowsMax<<
1870  ", NRowsInWhichStructVarArePresent = "<<
1871  NRowsInWhichStructVarArePresent<<", nRanges = "<<nRanges
1872  <<", nBounds = "<<nBounds<<endl;
1873 */
1874  int status= glp_main(
1875  nRows,nameRows,typeRows, // ROWS info
1876  NStructVar, NStructRowsMax, NRowsInWhichStructVarArePresent, // COLUMNS info
1877  StructVarName, NameRowsInWhichStructVarArePresent, // COLUMNS info
1878  Coefficients, // COLUMNS info
1879  ValueB, // RHS info
1880  nRanges, ValueRanges, NameRanges, // RANGES info
1881  nBounds, BoundValue, BoundStructVarName, TypeofBound // BOUNDS info
1882 // ,final_values, TIMEOUT
1883  ,final_values
1884  );
1885 
1886 
1887  if (status != 0) return -5; // fit failed
1888 
1889 //--------stampaggi
1890 /*
1891 printf("from main, final printout con routines chiamate direttamente -------------------------------\n");
1892 printf(" number of structural variables %d\n",NStructVar);
1893 int ica;
1894 for(ica=0;ica<NStructVar;ica++){
1895  printf("name of structural variable %s and its final value %g\n",
1896  StructVarName[ica], final_values[ica]);
1897 }
1898 printf("from main, end of final printout con routines chiamate direttamente -------------------------------\n");
1899 */
1900 //--------fine stampaggi
1901 
1902 
1903 
1904 //----------------------- fine funzioni chiamate direttamente
1905 
1906 
1907  m1_result=final_values[0];
1908  m2_result=final_values[1];
1909 // q1_result=final_values[2];
1910 // q2_result=final_values[3];
1911 
1912 
1913  *emme = m1_result-m2_result ;
1914 // taking into account the rotation + traslation that was performed and calculate emme and qu
1915 
1916  if(fabs(cose-*emme*sine)> 1.e-10) {
1917  *emme=((*emme)*cose+sine)/(cose-(*emme)*sine);
1918  return 1;
1919  } else { // in this case the equation is 0 = U in the Conformal plane --> x=0 in the XY plane.
1920  return -99;
1921  }
1922 
1923 
1924 
1925 }
1926 
1927 //----------end of function PndTrkGlpkFits::FitSZspace
1928 
1929 
1930 
1931 
1932 ClassImp(PndTrkGlpkFits)
1933 
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
i
Int_t i
Definition: run_full.C:25
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
n
int n
Definition: HypStatDecay.cxx:21
PndTrkGlpkFits::FitSZspace
Short_t FitSZspace(Short_t nSkewHitsinTrack, Double_t *S, Double_t *Z, Double_t *DriftRadius, Double_t *ErrorDriftRadius, Double_t FInot, Short_t NMAX, Double_t *emme, int IVOLTE)
Definition: PndTrkGlpkFits.cxx:986
PndTrkGlpkFits::FitHelixCylinder
Short_t FitHelixCylinder(Short_t nHitsinTrack, Double_t *Xconformal, Double_t *Yconformal, Double_t *DriftRadiusconformal, Double_t *ErrorDriftRadiusconformal, Double_t rotationangle, Double_t trajectory_vertex[2], Short_t NMAX, Double_t *m, Double_t *q, Double_t *pAlfa, Double_t *pBeta, Double_t *pGamma, bool *Type, int istampa, int IVOLTE)
Definition: PndTrkGlpkFits.cxx:14
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PI
#define PI
Definition: createSttGeometry.C:9
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
ClassImp
ClassImp(PndAnaContFact)
Z
double Z
Definition: anaLmdDigi.C:68
angle
Double_t angle
Definition: createROOTEnDrcGeo.C:13
status
int status[10]
Definition: f_Init.h:28
Generated on Wed Apr 3 2019 10:02:59 for FairRoot/PandaRoot by   doxygen 1.8.5