padsley / k600analyser

Code for the K600 analyser including plugin codes for silicon, clover and HAGAR data
1 stars 4 forks source link

Resolution online and offline #141

Closed padsley closed 8 years ago

padsley commented 8 years ago

Currently the resolution online and offline differ: the resolution online seems much better than in the offline analysis.

Why? How do we fix this?

The online version of 'raytrace' is shown below.

void raytrace(Double_t dd[],Int_t wire[],Double_t *_X,Double_t *_Th,Double_t *_chisq,Int_t wire_num,Double_t res[],Int_t *_flag){

   //printf("-----------------start raytrace-----------------\n");
   int i;
   Double_t x_tttv, sum_0=0.0, sum_x=0.0, sum_z=0.0, sum_xz=0.0, sum_x2=0.0;
   Int_t    wireID_first=0,wireID_last=0,wireID_min=0;
   Double_t wireID_ave=0,wirechisq=0;
   Double_t x_first, x_last;
   Double_t a1,b1,a2,b2,a,b;
   Double_t a_pre, b_pre, X_pre, Th_pre;
   Double_t X_v1=1.0, Th_v1=1.0,chisq_v1=0.0;
   Double_t X_v2=1.0, Th_v2=1.0,chisq_v2=0.0;
   Double_t drift_min_esti,drift_difference;
   Int_t flag1=0;
   Double_t tmpdd[20];
   Int_t tmpwire[20];

   wireID_first = 0;                // wireID_first an array index nr, NOT a real wire number
   wireID_last  = wire_num-1;       // wireID_last  an array index nr, NOT a real wire number

   // Label wire associated with minimum drift distance
   for(i=1;i<wire_num;i++){
     if(dd[wireID_min] > dd[i]) wireID_min = i;       //wireID_min an array index nr, NOT a real wire number
   }

   // Label events with more than 1 drift distance local minimum
   for(i=1;i<(wire_num-1);i++){
     if( (dd[i]>dd[i-1]) && (dd[i]>dd[i+1]) ) {
       flag1=1;    
       if(dd[i-1]==dd[wireID_first]){   // is one of the minima at start of wiregroup?
          flag1=2;
       }
       if(dd[i+1]==dd[wireID_last]){   // is one of the minima at the end of the wiregroup?
          flag1=3;
       }
     }
   }

   if(flag1!=0){        // if I know there are more than 2 min I try use only the middle wires
     for(i=1;i<(wire_num-1);i++){
       tmpdd[i-1]=dd[i];
       tmpwire[i-1]=wire[i];
       dd[i-1]=tmpdd[i-1];
       wire[i-1]=tmpwire[i-1];
     }
     wire_num=wire_num-2;
     wireID_last= wire_num-1;

     // TRY AGAIN: Label wire associated with minimum drift distance
     flag1=0;
     for(i=1;i<wire_num;i++){
       if(dd[wireID_min] > dd[i]) wireID_min = i;       //wireID_min an array index nr, NOT a real wire number
     }

     Int_t counter=1;
     // How to identify monotinically increasing   sets of dd
     for(i=0;i<(wire_num-1);i++){
       if( dd[i]>dd[i+1] ) {
     counter+=1;    
       }
     }    
     if(counter==wire_num) flag1=5;

     counter=1;
     // How to identify monotinically decreasing  sets of dd
     for(i=0;i<(wire_num-1);i++){
       if( dd[i]<dd[i+1] ) {
     counter+=1;    
       }
     }    
     if(counter==wire_num) flag1=6;

     // TRY AGAIN Label events with more than 1 drift distance local minimum
     for(i=1;i<(wire_num-1);i++){
       if( (dd[i]>dd[i-1]) && (dd[i]>dd[i+1]) ) {
     flag1=1;    
     if(dd[i-1]==dd[wireID_first]){   // is one of the minima at start of wiregroup?
       flag1=2;
     }
     if(dd[i+1]==dd[wireID_last]){   // is one of the minima at the end of the wiregroup?
       flag1=3;
     }
       }
     }
     if(wire_num<3) flag1=4;
   }   

   // Determine tentative ray: straigh line through first and last wires
   x_first = X_WIRESPACE*(wire[wireID_first]);  //pos in mm of first wire of group
   x_last  = X_WIRESPACE*(wire[wireID_last]);   //pos in mm of last wire of group

   a_pre = (dd[wireID_first]+dd[wireID_last])/(x_first-x_last);    // slope
   b_pre = dd[wireID_first]-a_pre*x_first;                        // y offset
   X_pre = (-1.)*b_pre/a_pre;
   Th_pre =(-1.)*57.29578*atan(a_pre);

   // Use the above information to determine sign of drift distance
   // + drift on K600 side, - drift downstream side
   for(i=0;i<wire_num;++i){     
     if( ( a_pre*(X_WIRESPACE*wire[i])+b_pre) < 0 ){
        dd[i] *= -1;
     }
   }

   //========================================================================
   //least square fit 1
   sum_0=0; sum_x=0; sum_z=0; sum_xz=0; sum_x2=0; 

   for(i=0;i<wire_num;++i){
     if(dd[i]!=dd[wireID_min]){         //ignore the wire with smallest drifttime. It causes inaccuracies
       x_tttv = X_WIRESPACE*(float(wire[i]));  // left -> high Ex  
       sum_0  += 1.0;
       sum_x  += x_tttv;
       sum_z  += dd[i];
       sum_xz += x_tttv*dd[i];
       sum_x2 += pow(x_tttv,2);
     }
   }
   a1 = (sum_x*sum_z-sum_0*sum_xz)/(pow(sum_x,2)-sum_0*sum_x2);
   b1 = (sum_x*sum_xz-sum_x2*sum_z)/(pow(sum_x,2)-sum_0*sum_x2);  
   X_v1  = (-1.)*b1/a1;
   Th_v1 = (-1.)*57.29578*atan(a1); // [deg.]

   for(i=0;i<wire_num;++i){
     if(dd[i]!=dd[wireID_min]){            //ignore the wire with smallest drifttime. It causes inaccuracies
       chisq_v1 += pow(dd[i]-(a1*X_WIRESPACE*(wire[i])+b1),2); 
       ++wirechisq;
     }
   }
   chisq_v1 /= wirechisq;

   //----- the values sent back to the main routine---------
   *_X = X_v1; 
   *_Th = Th_v1;             
   *_chisq = chisq_v1;
   a=a1;
   b=b1;
   *_flag=flag1;

   //----------------------------------------------

   for(i=0;i<7;i++){
     res[i]=-100.;     // initial value of res array 
   }

   if(wireID_first!=wireID_min && wireID_last!=wireID_min && wire_num>3){  
   // min drift distance may not be on edge of wire group, and only look at events >4 wires
     //res[0] is a relative position in drift cell where track cross wireplane
     res[0]=(*_X/4.0)-int(*_X/4.);           

     //res[1] same as XSYS calc: the diff in slope (since wireID_min is used res[1] not very accurate
     res[1]=(dd[wireID_min+1]-dd[wireID_min])/(wire[wireID_min+1]-wire[wireID_min]) 
       + ((dd[wireID_min-1]-dd[wireID_min])/(wire[wireID_min]-wire[wireID_min-1]));  

     //res[5] is difference in pos calculated at min driftt wire; drift_min_esti from RCNP analyzer sorter_user.c
     drift_min_esti=( (wireID_last-wireID_min)*dd[wireID_first] + (wireID_min-wireID_first)*dd[wireID_last] )
                    / (wireID_last - wireID_first);   
     res[5]=drift_min_esti- dd[wireID_min];

     //res[6] is close to the Bertozzi definition
     if((wireID_last-wireID_first+1)>4 && wire_num>4) {
       res[6]=(dd[wireID_first]-dd[wireID_first+1]) - (dd[wireID_last-1] - dd[wireID_last]);
     }
     else res[6]=-9.;

     //res[7] is drift_min_esti from RCNP analyzer sorter_user.c
     res[7]=drift_min_esti;

     // drift_min_esti from EVAL code?????true???
     drift_min_esti = a* (X_WIRESPACE*(wire[wireID_min])) + b;   
     res[8]=drift_min_esti;
   }

}
padsley commented 8 years ago

My plan is to add a 'raytrace_old' function to one of the analysers on arthur so that I can test resolution effects on that computer.

padsley commented 8 years ago

So, we think that (having looked a bit more into this) that the problem is actually not in the sorting but is rather in the beam behaviour not being consistent during the experiment. While CAKE was put in, the fields have shifted and screwed up the resolution. :(