/[PAMELA software]/DarthVader/ToFLevel2/src/ToFLevel2.cpp
ViewVC logotype

Diff of /DarthVader/ToFLevel2/src/ToFLevel2.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.18 by mocchiut, Mon Nov 26 08:01:17 2007 UTC revision 1.25 by pamelats, Wed Apr 15 10:05:22 2009 UTC
# Line 1  Line 1 
1  /**  /**
2   * \file ToFLevel2.cpp   * \file ToFLevel2.cpp
3   * \author Gianfranca DeRosa, Wolfgang Menn   * \author Gianfranca DeRosa, Wolfgang Menn
4     *
5     * WM dec 2008: Description of "GetdEdx" changed
6     * WM dec 2008: "GetdEdxPaddle" modified: Now includes saturation limit
7     *              PMTs higher than the saturation limit are not used for dEdx
8     * WM apr 2009: bug found by Nicola in method "GetPaddlePlane"
9   */   */
10    
 #include <TObject.h>  
11  #include <ToFLevel2.h>  #include <ToFLevel2.h>
 #include <iostream>  
12  using namespace std;  using namespace std;
13  ClassImp(ToFPMT);  ClassImp(ToFPMT);
14  ClassImp(ToFTrkVar);  ClassImp(ToFTrkVar);
# Line 216  Int_t ToFLevel2::GetNHitPaddles(Int_t pl Line 219  Int_t ToFLevel2::GetNHitPaddles(Int_t pl
219      return npad;      return npad;
220  };  };
221    
222    //wm Nov 08
223  //gf Apr 07  //gf Apr 07
224  /**  /**
225   * Method to get the mean dEdx from a given ToF plane. This current version   * Method to get the mean dEdx from a ToF layer - ATTENTION:
226   * is just summing up all PMT signals, which will not give proper results,   * It will sum up the dEdx of all the paddles, but since by definition
227   *  and needs a revision.   * only the paddle hitted by the track gets a dEdx value and the other
228     * paddles are set to zero, the output is just the dEdx of the hitted
229     * paddle in each layer!
230     * The "adcfl" option is not very useful (an artificial dEdx is per
231     * definition= 1 mip and not a real measurement), anyway left in the code
232   * @param notrack Track Number   * @param notrack Track Number
233   * @param plane Plane index (0,1,2,3,4,5)   * @param plane Plane index (0,1,2,3,4,5)
234   * @param adcflag in the plane (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )   * @param adcflag in the plane (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )
# Line 381  void ToFLevel2::GetPMTIndex(Int_t ind, I Line 388  void ToFLevel2::GetPMTIndex(Int_t ind, I
388    
389    
390    
391    //  wm Nov 08 revision - saturation values included
392  /// gf Apr 07  /// gf Apr 07
   
393  /**  /**
394   * Method to get the dEdx from a given ToF paddle.   * Method to get the dEdx from a given ToF paddle.
395     * If two PMTs are good, the mean dEdx of both PMTs is taken, otherwise
396     * just the dEdx of the "good" PMT. If both PMTs are above saturation => dEdx=1000
397   * @param notrack Track Number   * @param notrack Track Number
398   * @param Paddle index (0,1,...,23).   * @param Paddle index (0,1,...,23).
399   * @param adcflag in the paddle (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )   * @param adcflag in the paddle (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )
# Line 394  void ToFLevel2::GetPMTIndex(Int_t ind, I Line 402  void ToFLevel2::GetPMTIndex(Int_t ind, I
402   */   */
403  void ToFLevel2::GetdEdxPaddle(Int_t notrack, Int_t paddleid, Int_t adcfl, Float_t &PadEdx, Int_t &SatWarning){  void ToFLevel2::GetdEdxPaddle(Int_t notrack, Int_t paddleid, Int_t adcfl, Float_t &PadEdx, Int_t &SatWarning){
404    
405    /*
406    Float_t  PMTsat[48] = {
407    3162.14, 3165.48, 3153.85, 3085.73, 3089.65, 3107.64, 3097.52, 3078.37,
408    3130.05, 3087.07, 3112.22, 3102.92, 3080.58, 3092.55, 3087.94, 3125.03,
409    3094.09, 3143.16, 3125.51, 3181.27, 3092.09, 3124.98, 3069.3, 3095.53,
410    3097.11, 3133.53, 3114.73, 3113.01, 3091.19, 3097.99, 3033.84, 3134.98,
411    3081.37, 3111.04, 3066.77, 3108.17, 3133, 3111.06, 3052.52, 3140.66,
412    3106.33, 3094.85, 3150.85, 3118.8, 3096.24, 3118.47,3111.36, 3117.11 } ;
413    */
414    
415    // new values from Napoli dec 2008
416    Float_t  PMTsat[48] = {
417    3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
418    3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
419    3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
420    3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
421    3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
422    3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };
423    
424    for (Int_t i=0; i<48;i++) PMTsat[i] = PMTsat[i] - 5.;  // safety margin
425    
426    
427    PadEdx = 0.;    PadEdx = 0.;
428    SatWarning = 1000;  //  SatWarning = 1000;
429      SatWarning = 0;   // 0=good, increase for each bad PMT
430    
431    Float_t dEdx[48] = {0};    Float_t dEdx[48] = {0};
432    Int_t pmt_id = -1;    Int_t pmt_id = -1;
# Line 427  void ToFLevel2::GetdEdxPaddle(Int_t notr Line 458  void ToFLevel2::GetdEdxPaddle(Int_t notr
458        adcraw[pmtright] = pmt->adc;        adcraw[pmtright] = pmt->adc;
459      }      }
460    }    }
461    
462        
463    for (Int_t i=0; i<trk->npmtadc; i++){    for (Int_t i=0; i<trk->npmtadc; i++){
464    
# Line 439  void ToFLevel2::GetdEdxPaddle(Int_t notr Line 471  void ToFLevel2::GetdEdxPaddle(Int_t notr
471      }      }
472    }    }
473    
474    if( adcraw[pmtleft] >3000 || adcraw[pmtright] >3000)SatWarning=1;  
475        //  if( adcraw[pmtleft] >3000 || adcraw[pmtright] >3000)SatWarning=1;  //old version
476    if(dEdx[pmtleft]!=0 && dEdx[pmtright]!=0){  
477      PadEdx = (dEdx[pmtleft]+dEdx[pmtright])*0.5;  // Increase SatWarning Counter for each PMT>Sat
478    }    if( adcraw[pmtleft] > PMTsat[pmtleft])SatWarning++;  
479    if(dEdx[pmtleft]==0 && dEdx[pmtright]!=0){    if( adcraw[pmtright] > PMTsat[pmtright])SatWarning++;
480      PadEdx = dEdx[pmtright];  
481    }  // if ADC  > sat set dEdx=1000
482    if(dEdx[pmtleft]!=0 && dEdx[pmtright]==0){    if( adcraw[pmtleft] > PMTsat[pmtleft]) dEdx[pmtleft] = 1000.;
483      PadEdx = dEdx[pmtleft];    if( adcraw[pmtright] > PMTsat[pmtright]) dEdx[pmtright] = 1000. ;
484    }  
485    // if two PMT are good, take mean dEdx, otherwise only the good dEdx
486      if(dEdx[pmtleft]<1000 && dEdx[pmtright]<1000) PadEdx = (dEdx[pmtleft]+dEdx[pmtright])*0.5;
487      if(dEdx[pmtleft]==1000 && dEdx[pmtright]<1000) PadEdx = dEdx[pmtright];  
488      if(dEdx[pmtleft]<1000 && dEdx[pmtright]==1000) PadEdx = dEdx[pmtleft];
489        
   return;  
490  };  };
491  //  //
492    
# Line 508  TString ToFLevel2::GetPMTName(Int_t ind) Line 543  TString ToFLevel2::GetPMTName(Int_t ind)
543        
544  };  };
545    
546    // wm jun 08
 // gf Apr 07  
547  Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane){  Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane){
548    return GetPaddleIdOfTrack(xtr ,ytr ,plane, 0.4);
549    }
550    
551    // gf Apr 07
552    Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane, Float_t margin){
553      
554    Double_t xt,yt,xl,xh,yl,yh;    Double_t xt,yt,xl,xh,yl,yh;
555        
556    Float_t tof11_x[8] = {-17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85};    Float_t tof11_x[8] = {-17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85};
# Line 541  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 580  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
580      yh =  33.0/2. ;      yh =  33.0/2. ;
581      if ((yt>yl)&&(yt<yh)) {      if ((yt>yl)&&(yt<yh)) {
582        for (Int_t i1=0; i1<8;i1++){        for (Int_t i1=0; i1<8;i1++){
583          xl = tof11_x[i1] - (5.1-0.4)/2. ;          xl = tof11_x[i1] - (5.1-margin)/2. ;
584          xh = tof11_x[i1] + (5.1-0.4)/2. ;          xh = tof11_x[i1] + (5.1-margin)/2. ;
585          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
586        }        }
587      }      }
# Line 559  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 598  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
598            
599      if ((xt>xl)&&(xt<xh)) {      if ((xt>xl)&&(xt<xh)) {
600        for (Int_t i1=0; i1<6;i1++){        for (Int_t i1=0; i1<6;i1++){
601          yl = tof12_y[i1] - (5.5-0.4)/2. ;          yl = tof12_y[i1] - (5.5-margin)/2. ;
602          yh = tof12_y[i1] + (5.5-0.4)/2. ;          yh = tof12_y[i1] + (5.5-margin)/2. ;
603          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;
604        }        }
605      }      }
# Line 577  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 616  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
616            
617      if ((xt>xl)&&(xt<xh)) {      if ((xt>xl)&&(xt<xh)) {
618        for (Int_t i1=0; i1<2;i1++){        for (Int_t i1=0; i1<2;i1++){
619          yl = tof21_y[i1] - (7.5-0.4)/2. ;          yl = tof21_y[i1] - (7.5-margin)/2. ;
620          yh = tof21_y[i1] + (7.5-0.4)/2. ;          yh = tof21_y[i1] + (7.5-margin)/2. ;
621          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;
622        }        }
623      }      }
# Line 594  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 633  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
633            
634      if ((yt>yl)&&(yt<yh)) {      if ((yt>yl)&&(yt<yh)) {
635        for (Int_t i1=0; i1<2;i1++){        for (Int_t i1=0; i1<2;i1++){
636          xl = tof22_x[i1] - (9.0-0.4)/2. ;          xl = tof22_x[i1] - (9.0-margin)/2. ;
637          xh = tof22_x[i1] + (9.0-0.4)/2. ;          xh = tof22_x[i1] + (9.0-margin)/2. ;
638          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
639        }        }
640      }      }
# Line 611  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 650  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
650            
651      if ((yt>yl)&&(yt<yh)) {      if ((yt>yl)&&(yt<yh)) {
652        for (Int_t i1=0; i1<3;i1++){        for (Int_t i1=0; i1<3;i1++){
653          xl = tof31_x[i1] - (6.0-0.4)/2. ;          xl = tof31_x[i1] - (6.0-margin)/2. ;
654          xh = tof31_x[i1] + (6.0-0.4)/2. ;          xh = tof31_x[i1] + (6.0-margin)/2. ;
655          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
656        }        }
657      }      }
# Line 628  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 667  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
667            
668      if ((xt>xl)&&(xt<xh)) {      if ((xt>xl)&&(xt<xh)) {
669        for (Int_t i1=0; i1<3;i1++){        for (Int_t i1=0; i1<3;i1++){
670          yl = tof32_y[i1] - (5.0-0.4)/2. ;          yl = tof32_y[i1] - (5.0-margin)/2. ;
671          yh = tof32_y[i1] + (5.0-0.4)/2. ;          yh = tof32_y[i1] + (5.0-margin)/2. ;
672          if ((yt>yl)&&(yt<yh)) paddleidoftrack=i1;          if ((yt>yl)&&(yt<yh)) paddleidoftrack=i1;
673        }        }
674      }      }
# Line 696  void ToFLevel2::GetPMTPaddle(Int_t pmt_i Line 735  void ToFLevel2::GetPMTPaddle(Int_t pmt_i
735  // gf Apr 07  // gf Apr 07
736    
737  void ToFLevel2::GetPaddlePMT(Int_t paddle, Int_t &pmtleft, Int_t &pmtright){  void ToFLevel2::GetPaddlePMT(Int_t paddle, Int_t &pmtleft, Int_t &pmtright){
738      pmtleft=paddle*2;
739    if(paddle==0){    pmtright= pmtleft+1;  
     pmtleft=0;  
     pmtright=1;  
   }  
   
   if(paddle==1){  
     pmtleft=2;  
     pmtright=3;  
   }  
   
   if(paddle==2){  
     pmtleft=4;  
     pmtright=5;  
   }  
   
   if(paddle==3){  
     pmtleft=6;  
     pmtright=7;  
   }  
   
   if(paddle==4){  
     pmtleft=8;  
     pmtright=9;  
   }  
   
   if(paddle==5){  
     pmtleft=10;  
     pmtright=11;  
   }  
   
   if(paddle==6){  
     pmtleft=12;  
     pmtright=13;  
   }  
   
   if(paddle==7){  
     pmtleft=14;  
     pmtright=15;  
   }  
   
   if(paddle==8){  
     pmtleft=16;  
     pmtright=17;  
   }  
   
   if(paddle==9){  
     pmtleft=18;  
     pmtright=19;  
   }  
   
   if(paddle==10){  
     pmtleft=20;  
     pmtright=21;  
   }  
   
   if(paddle==11){  
     pmtleft=22;  
     pmtright=23;  
   }  
   
   if(paddle==12){  
     pmtleft=24;  
     pmtright=25;  
   }  
   
   if(paddle==13){  
     pmtleft=26;  
     pmtright=27;  
   }  
   
   if(paddle==14){  
     pmtleft=28;  
     pmtright=29;  
   }  
   
   if(paddle==15){  
     pmtleft=30;  
     pmtright=31;  
   }  
   
   if(paddle==16){  
     pmtleft=32;  
     pmtright=33;  
   }  
   
   if(paddle==17){  
     pmtleft=34;  
     pmtright=35;  
   }  
   
   if(paddle==18){  
     pmtleft=36;  
     pmtright=37;  
   }  
   
   if(paddle==19){  
     pmtleft=38;  
     pmtright=39;  
   }  
   
   if(paddle==20){  
     pmtleft=40;  
     pmtright=41;  
   }  
   
   if(paddle==21){  
     pmtleft=42;  
     pmtright=43;  
   }  
   
   if(paddle==22){  
     pmtleft=44;  
     pmtright=45;  
   }  
   
   if(paddle==23){  
     pmtleft=46;  
     pmtright=47;  
   }  
     
740    return;    return;
741  }  }
742    
# Line 930  void ToFLevel2::GetPaddleGeometry(Int_t Line 850  void ToFLevel2::GetPaddleGeometry(Int_t
850   */   */
851  Int_t ToFLevel2::GetPaddleid(Int_t plane, Int_t paddle)  Int_t ToFLevel2::GetPaddleid(Int_t plane, Int_t paddle)
852  {  {
   
853    Int_t padid=-1;    Int_t padid=-1;
854    Int_t pads11=8;    Int_t pads[6]={8,6,2,2,3,3};
   Int_t pads12=6;  
   Int_t pads21=2;  
   Int_t pads22=2;  
   Int_t pads31=3;  
   //  Int_t pads32=3;  
   
   
   if(plane == 0){  
     padid=paddle;  
   }  
   
   if(plane == 1){  
     padid=pads11+paddle;  
   }  
   
   if(plane == 2){  
     padid=pads11+pads12+paddle;  
   }  
   
   if(plane == 3){  
     padid=pads11+pads12+pads21+paddle;  
   }  
   
   if(plane == 4){  
     padid=pads11+pads12+pads21+pads22+paddle;  
   }  
855    
856    if(plane == 5){    int somma=0;
857      padid=pads11+pads12+pads21+pads22+pads31+paddle;    int np=plane;
858      for(Int_t j=0; j<np; j++){
859        somma+=pads[j];
860    }    }
861      padid=paddle+somma;
862    return padid;    return padid;
863    
864  }  }
# Line 994  void ToFLevel2::GetPaddlePlane(Int_t pad Line 889  void ToFLevel2::GetPaddlePlane(Int_t pad
889      return;      return;
890    }    }
891    
892    if(7<pad<14){    if((7<pad)&&(pad<14)){
893      plane=1;      plane=1;
894      paddle=pad-pads11;      paddle=pad-pads11;
895      return;      return;
896    }    }
897        
898    if(13<pad<16){    if((13<pad)&&(pad<16)){
899      plane=2;      plane=2;
900      paddle=pad-pads11-pads12;      paddle=pad-pads11-pads12;
901      return;      return;
902    }    }
903    
904    if(15<pad<18){    if((15<pad)&&(pad<18)){
905      plane=3;      plane=3;
906      paddle=pad-pads11-pads12-pads21;      paddle=pad-pads11-pads12-pads21;
907      return;      return;
908    }    }
909    
910    if(17<pad<21){    if((17<pad)&&(pad<21)){
911      plane=4;      plane=4;
912      paddle=pad-pads11-pads12-pads21-pads22;      paddle=pad-pads11-pads12-pads21-pads22;
913      return;      return;
914    }    }
915    
916    if(20<pad<24){    if((20<pad)&&(pad<24)){
917      plane=5;      plane=5;
918      paddle=pad-pads11-pads12-pads21-pads22-pads31;      paddle=pad-pads11-pads12-pads21-pads22-pads31;
919      return;      return;
# Line 1049  Int_t ToFLevel2::GetNPaddle(Int_t plane) Line 944  Int_t ToFLevel2::GetNPaddle(Int_t plane)
944    
945  }  }
946    
 ////////////////////////////////////////////////////  
947    
948    
949    /// wm feb 08
950    
951    /**
952     * Method to calculate Beta from the 12 single measurements
953     * we check the individual weights for artificial TDC values, then calculate
954     * am mean beta for the first time. In a second step we loop again through
955     * the single measurements, checking for the residual from the mean
956     * The cut on the residual reject measurements > "x"-sigma. A chi2 value is
957     * calculated, furthermore a "quality" value by adding the weights which
958     * are finally used. If all measurements are taken, "quality" will be = 22.47.
959     * A chi2 cut around 3-4 and a quality-cut > 20 is needed for clean beta
960     * measurements like antiprotons etc.
961     * The Level2 output is derived in the fortran routines using: 10.,10.,20.
962     * @param notrack Track Number
963     * @param cut on residual: difference between single measurement and mean
964     * @param cut on "quality"
965     * @param cut on chi2
966     */
967    
968    Float_t ToFLevel2::CalcBeta(Int_t notrack, Float_t resmax, Float_t qualitycut, Float_t chi2cut){
969    
970    //  cout<<" in CalcBeta "<<resmax<<" "<<chi2cut<<" "<<qualitycut<<endl;
971    
972      Float_t bxx = 100.;
973      //
974      ToFTrkVar *trk = GetToFTrkVar(notrack);
975      if(!trk) return 0; //ELENA
976    
977    
978      Float_t chi2,xhelp,beta_mean;
979      Float_t w_i[12],quality,sw,sxw,res,betachi,beta_mean_inv;
980      Float_t b[12],tdcfl;
981      Int_t  pmt_id,pmt_plane;
982    
983      for (Int_t i=0; i<12; i++){
984        b[i] = trk->beta[i];
985                                  }
986          
987    
988    //========================================================================
989    //---  Find out ToF layers with artificial TDC values & fill vector    ---
990    //========================================================================
991    
992    Float_t  w_il[6];
993    
994         for (Int_t jj=0; jj<6;jj++) {
995             w_il[jj] = 1000.;
996                                     }
997    
998    
999      for (Int_t i=0; i<trk->npmttdc; i++){
1000        //
1001        pmt_id = (trk->pmttdc).At(i);
1002        pmt_plane = GetPlaneIndex(pmt_id);
1003        tdcfl = (trk->tdcflag).At(i);
1004        if (w_il[pmt_plane] != 1.) w_il[pmt_plane] = tdcfl; //tdcflag
1005                                         };
1006      
1007    //========================================================================
1008    //---  Set weights for the 12 measurements using information for top and bottom:
1009    //---  if no measurements: weight = set to very high value=> not used
1010    //---  top or bottom artificial: weight*sqrt(2)
1011    //---  top and bottom artificial: weight*sqrt(2)*sqrt(2)
1012    //========================================================================
1013    
1014    Int_t itop[12] = {0,0,1,1,2,2,3,3,0,0,1,1};
1015    Int_t ibot[12] = {4,5,4,5,4,5,4,5,2,3,2,3};
1016    
1017         xhelp= 1E09;
1018      
1019         for (Int_t jj=0; jj<12;jj++) {
1020         if (jj<4)           xhelp = 0.11;    // S1-S3
1021         if ((jj>3)&&(jj<8)) xhelp = 0.18;    // S2-S3
1022         if (jj>7)           xhelp = 0.28;    // S1-S2
1023         if ((w_il[itop[jj]] == 1000.) && (w_il[ibot[jj]] == 1000.)) xhelp = 1E09;
1024         if ((w_il[itop[jj]] == 1) || (w_il[ibot[jj]] == 1.)) xhelp = xhelp*1.414 ;
1025         if ((w_il[itop[jj]] == 1) && (w_il[ibot[jj]] == 1.)) xhelp = xhelp*2. ;
1026    
1027         w_i[jj] = 1./xhelp;
1028                                      }
1029    
1030    
1031    //========================================================================
1032    //--- Calculate mean beta for the first time -----------------------------
1033    //--- We are using "1/beta" since its error is gaussian ------------------
1034    //========================================================================
1035    
1036          Int_t icount=0;
1037          sw=0.;
1038          sxw=0.;
1039          beta_mean=100.;
1040    
1041              for (Int_t jj=0; jj<12;jj++){
1042            if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.))
1043             {
1044                icount= icount+1;
1045                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1046                sw =sw + w_i[jj]*w_i[jj] ;
1047    
1048             }
1049             }
1050    
1051          if (icount>0) beta_mean=1./(sxw/sw);
1052          beta_mean_inv = 1./beta_mean;
1053    
1054    //========================================================================
1055    //--- Calculate beta for the second time, use residuals of the single
1056    //--- measurements to get a chi2 value
1057    //========================================================================
1058    
1059          icount=0;
1060          sw=0.;
1061          sxw=0.;
1062          betachi = 100.;
1063          chi2 = 0.;
1064          quality=0.;
1065    
1066    
1067              for (Int_t jj=0; jj<12;jj++){
1068           if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.)&&(w_i[jj]>0.01)) {
1069                res = beta_mean_inv - (1./b[jj]) ;
1070                if (fabs(res*w_i[jj])<resmax)          {;
1071                chi2 = chi2 + pow((res*w_i[jj]),2) ;
1072                icount= icount+1;
1073                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1074                sw =sw + w_i[jj]*w_i[jj] ;
1075                                                   }
1076                                                                            }
1077                                          }
1078          quality = sqrt(sw) ;
1079    
1080          if (icount==0) chi2 = 1000.;
1081          if (icount>0) chi2 = chi2/(icount) ;
1082          if (icount>0) betachi=1./(sxw/sw);
1083    
1084       bxx = 100.;
1085       if ((chi2 < chi2cut)&&(quality>qualitycut)) bxx = betachi;
1086      //
1087      return(bxx);
1088    };
1089    
1090    
1091    ////////////////////////////////////////////////////
1092    ////////////////////////////////////////////////////
1093    
1094    
1095  /**  /**
1096   * Fills a struct cToFLevel2 with values from a ToFLevel2 object (to put data into a F77 common).   * Fills a struct cToFLevel2 with values from a ToFLevel2 object (to put data into a F77 common).
# Line 1099  void ToFLevel2::GetLevel2Struct(cToFLeve Line 1138  void ToFLevel2::GetLevel2Struct(cToFLeve
1138        }        }
1139    } //ELENA    } //ELENA
1140  }  }
1141    
1142    
1143    //
1144    // Reprocessing tool // Emiliano 08/04/07
1145    //
1146    Int_t ToFLevel2::Process(TrkLevel2 *trk, TrigLevel2 *trg, GL_RUN *run, OrbitalInfo *orb, Bool_t force){
1147      //
1148      // Copiare qui qualcosa di simile a calonuclei per evitare di riprocessare sempre tutto
1149      //
1150    
1151    
1152    
1153    
1154      //
1155      // structures to communicate with F77
1156      //
1157      extern struct ToFInput  tofinput_;
1158      extern struct ToFOutput tofoutput_;
1159      //
1160      // DB connection
1161      //
1162      TString host;
1163      TString user;
1164      TString psw;
1165      const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
1166      const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
1167      const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
1168      if ( !pamdbhost ) pamdbhost = "";
1169      if ( !pamdbuser ) pamdbuser = "";
1170      if ( !pamdbpsw ) pamdbpsw = "";
1171      if ( strcmp(pamdbhost,"") ) host = pamdbhost;
1172      if ( strcmp(pamdbuser,"") ) user = pamdbuser;
1173      if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
1174      //
1175      //
1176      TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
1177      if ( !dbc->IsConnected() ) return 1;
1178      stringstream myquery;
1179      myquery.str("");
1180      myquery << "SET time_zone='+0:00'";
1181      dbc->Query(myquery.str().c_str());
1182      GL_PARAM *glparam = new GL_PARAM();
1183      glparam->Query_GL_PARAM(1,1,dbc); // parameters stored in DB in GL_PRAM table
1184      trk->LoadField(glparam->PATH+glparam->NAME);
1185      //
1186      Bool_t defcal = true;
1187      Int_t error=glparam->Query_GL_PARAM(run->RUNHEADER_TIME,201,dbc); // parameters stored in DB in GL_PRAM table
1188      if ( error<0 ) {
1189        return(1);
1190      };
1191      printf(" Reading ToF parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
1192      if ( (UInt_t)glparam->TO_TIME != (UInt_t)4294967295UL ) defcal = false;
1193      //
1194      Int_t nlen = (Int_t)(glparam->PATH+glparam->NAME).Length();
1195      rdtofcal((char *)(glparam->PATH+glparam->NAME).Data(),&nlen);
1196      //
1197      Int_t adc[4][12];
1198      Int_t tdc[4][12];
1199      Float_t tdcc[4][12];
1200      //
1201      // process tof data
1202      //
1203      for (Int_t hh=0; hh<12;hh++){
1204        for (Int_t kk=0; kk<4;kk++){
1205               adc[kk][hh] = 4095;
1206               tdc[kk][hh] = 4095;
1207               tdcc[kk][hh] = 4095.;
1208               tofinput_.adc[hh][kk] = 4095;
1209               tofinput_.tdc[hh][kk] = 4095;
1210        };
1211      };
1212      Int_t ntrkentry = 0;
1213      Int_t npmtentry = 0;
1214      Int_t gg = 0;
1215      Int_t hh = 0;
1216      Int_t adcf[48];
1217      memset(adcf, 0, 48*sizeof(Int_t));
1218      Int_t tdcf[48];
1219      memset(tdcf, 0, 48*sizeof(Int_t));
1220      for (Int_t pm=0; pm < this->ntrk() ; pm++){
1221         ToFTrkVar *ttf = this->GetToFTrkVar(pm);
1222         for ( Int_t nc=0; nc < ttf->npmttdc; nc++){
1223                if ( (ttf->tdcflag).At(nc) != 0 ) tdcf[(ttf->pmttdc).At(nc)] = 1;
1224         };
1225         for ( Int_t nc=0; nc < ttf->npmtadc; nc++){
1226                if ( (ttf->adcflag).At(nc) != 0 ) adcf[(ttf->pmtadc).At(nc)] = 1;
1227         };
1228      };
1229      //
1230      for (Int_t pm=0; pm < this->npmt() ; pm++){
1231         ToFPMT *pmt = this->GetToFPMT(pm);
1232         this->GetPMTIndex(pmt->pmt_id, gg, hh);
1233         if ( adcf[pmt->pmt_id] == 0 ){
1234                 tofinput_.adc[gg][hh] = (int)pmt->adc;
1235                 adc[hh][gg] = (int)pmt->adc;
1236         };
1237         if ( tdcf[pmt->pmt_id] == 0 ){
1238                 tofinput_.tdc[gg][hh] = (int)pmt->tdc;
1239                 tdc[hh][gg] = (int)pmt->tdc;
1240         };
1241         tdcc[hh][gg] = (float)pmt->tdc_tw;
1242         // Int_t pppid = this->GetPMTid(hh,gg);
1243         //      printf(" pm %i pmt_id %i pppid %i hh %i gg %i tdcc %f tdc %f adc %f \n",pm,pmt->pmt_id,pppid,hh,gg,pmt->tdc_tw,pmt->tdc,pmt->adc);
1244      };
1245      //
1246      Int_t unpackError = this->unpackError;
1247      //
1248      for (Int_t hh=0; hh<5;hh++){
1249         tofinput_.patterntrig[hh]=trg->patterntrig[hh];
1250      };
1251      //
1252      this->Clear();
1253      //
1254          Int_t pmt_id = 0;
1255          ToFPMT *t_pmt = new ToFPMT();
1256          if(!(this->PMT)) this->PMT = new TClonesArray("ToFPMT",12); //ELENA
1257          TClonesArray &tpmt = *this->PMT;
1258          ToFTrkVar *t_tof = new ToFTrkVar();
1259          if(!(this->ToFTrk)) this->ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
1260          TClonesArray &t = *this->ToFTrk;
1261          //
1262          //
1263          // Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract common and track-related  variables.
1264          //
1265          npmtentry = 0;
1266          //
1267          ntrkentry = 0;
1268          //
1269          // Calculate tracks informations from ToF alone
1270          //
1271          tofl2com();
1272          //
1273          memcpy(this->tof_j_flag,tofoutput_.tof_j_flag,6*sizeof(Int_t));
1274          //
1275          t_tof->trkseqno = -1;
1276          //
1277          // and now we must copy from the output structure to the level2 class:
1278          //
1279          t_tof->npmttdc = 0;
1280          //
1281          for (Int_t hh=0; hh<12;hh++){
1282            for (Int_t kk=0; kk<4;kk++){
1283              if ( tofoutput_.tofmask[hh][kk] != 0 ){
1284                pmt_id = this->GetPMTid(kk,hh);
1285                t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
1286                t_tof->tdcflag.AddAt(tofoutput_.tdcflagtof[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
1287                t_tof->npmttdc++;
1288              };
1289            };
1290          };
1291          for (Int_t kk=0; kk<13;kk++){
1292            t_tof->beta[kk] = tofoutput_.betatof_a[kk];
1293          }
1294          //
1295          t_tof->npmtadc = 0;
1296          for (Int_t hh=0; hh<12;hh++){
1297            for (Int_t kk=0; kk<4;kk++){
1298              if ( tofoutput_.adctof_c[hh][kk] < 1000 ){
1299                t_tof->dedx.AddAt(tofoutput_.adctof_c[hh][kk],t_tof->npmtadc);
1300                pmt_id = this->GetPMTid(kk,hh);
1301                t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
1302                t_tof->adcflag.AddAt(tofoutput_.adcflagtof[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
1303                t_tof->npmtadc++;
1304              };
1305            };
1306          };
1307          //
1308          memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
1309          memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
1310          memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
1311          memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
1312          //
1313          new(t[ntrkentry]) ToFTrkVar(*t_tof);
1314          ntrkentry++;
1315          t_tof->Clear();
1316          //
1317          //
1318          //
1319          t_pmt->Clear();
1320          //
1321          for (Int_t hh=0; hh<12;hh++){
1322            for (Int_t kk=0; kk<4;kk++){
1323             // new WM
1324              if ( tofoutput_.tdc_c[hh][kk] < 4095 || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){
1325    //          if ( tdcc[kk][hh] < 4095. || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){
1326                //
1327                t_pmt->pmt_id = this->GetPMTid(kk,hh);
1328                t_pmt->tdc_tw = tofoutput_.tdc_c[hh][kk];
1329                t_pmt->adc = (Float_t)adc[kk][hh];
1330                t_pmt->tdc = (Float_t)tdc[kk][hh];
1331                //
1332                new(tpmt[npmtentry]) ToFPMT(*t_pmt);
1333                npmtentry++;
1334                t_pmt->Clear();
1335              };
1336            };
1337          };
1338          //
1339          // Calculate track-related variables
1340          //
1341          if ( trk->ntrk() > 0 ){
1342            //
1343            // We have at least one track
1344            //
1345            //
1346            // Run over tracks
1347            //
1348            for(Int_t nt=0; nt < trk->ntrk(); nt++){
1349              //
1350              TrkTrack *ptt = trk->GetStoredTrack(nt);
1351              //
1352              // Copy the alpha vector in the input structure
1353              //
1354              for (Int_t e = 0; e < 5 ; e++){
1355                tofinput_.al_pp[e] = ptt->al[e];
1356              };
1357              //
1358              // Get tracker related variables for this track
1359              //
1360              toftrk();
1361              //
1362              // Copy values in the class from the structure (we need to use a temporary class to store variables).
1363              //
1364              t_tof->npmttdc = 0;
1365              for (Int_t hh=0; hh<12;hh++){
1366                for (Int_t kk=0; kk<4;kk++){
1367                  if ( tofoutput_.tofmask[hh][kk] != 0 ){
1368                    pmt_id = this->GetPMTid(kk,hh);
1369                    t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
1370                    t_tof->tdcflag.AddAt(tofoutput_.tdcflag[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
1371                    t_tof->npmttdc++;
1372                  };
1373                };
1374              };
1375              for (Int_t kk=0; kk<13;kk++){
1376                t_tof->beta[kk] = tofoutput_.beta_a[kk];
1377              };
1378              //
1379              t_tof->npmtadc = 0;
1380              for (Int_t hh=0; hh<12;hh++){
1381                for (Int_t kk=0; kk<4;kk++){
1382                  if ( tofoutput_.adc_c[hh][kk] < 1000 ){
1383                    t_tof->dedx.AddAt(tofoutput_.adc_c[hh][kk],t_tof->npmtadc);
1384                    pmt_id = this->GetPMTid(kk,hh);
1385                    t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
1386                    t_tof->adcflag.AddAt(tofoutput_.adcflag[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
1387                    t_tof->npmtadc++;
1388                  };
1389                };
1390              };
1391              //
1392              memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
1393              memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
1394              memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
1395              memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
1396              //
1397              // Store the tracker track number in order to be sure to have shyncronized data during analysis
1398              //
1399              t_tof->trkseqno = nt;
1400              //
1401              // create a new object for this event with track-related variables
1402              //
1403              new(t[ntrkentry]) ToFTrkVar(*t_tof);
1404              ntrkentry++;
1405              t_tof->Clear();
1406              //
1407            }; // loop on all the tracks
1408          //
1409          this->unpackError = unpackError;
1410          if ( defcal ){
1411            this->default_calib = 1;
1412          } else {
1413            this->default_calib = 0;
1414          };
1415     };
1416    
1417    
1418    
1419      return(0);
1420    }
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.23