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