/[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.26 by carbone, Fri Nov 20 11:05:21 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(ToFdEdx);
15    ClassImp(ToFGeom);
16  ClassImp(ToFTrkVar);  ClassImp(ToFTrkVar);
17  ClassImp(ToFLevel2);  ClassImp(ToFLevel2);
18    
# Line 25  ToFPMT::ToFPMT(const ToFPMT &t){ Line 30  ToFPMT::ToFPMT(const ToFPMT &t){
30    tdc = t.tdc;    tdc = t.tdc;
31  }  }
32    
33  void ToFPMT::Clear(){  void ToFPMT::Clear(Option_t *t){
34    pmt_id = 0;    pmt_id = 0;
35    adc = 0.;    adc = 0.;
36    tdc_tw = 0.;    tdc_tw = 0.;
# Line 53  ToFTrkVar::ToFTrkVar() { Line 58  ToFTrkVar::ToFTrkVar() {
58    //    //
59  };  };
60    
61  void ToFTrkVar::Clear() {  void ToFTrkVar::Clear(Option_t *t) {
62    trkseqno = 0;    trkseqno = 0;
63    npmttdc = 0;    npmttdc = 0;
64    npmtadc = 0;    npmtadc = 0;
# Line 107  void ToFLevel2::Set(){//ELENA Line 112  void ToFLevel2::Set(){//ELENA
112      if(!ToFTrk)ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA      if(!ToFTrk)ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
113  }//ELENA  }//ELENA
114    
115  void ToFLevel2::Clear(){  void ToFLevel2::Clear(Option_t *t){
116    //    //
117    if(ToFTrk)ToFTrk->Delete(); //ELENA    if(ToFTrk)ToFTrk->Delete(); //ELENA
118    if(PMT)PMT->Delete(); //ELENA    if(PMT)PMT->Delete(); //ELENA
# Line 116  void ToFLevel2::Clear(){ Line 121  void ToFLevel2::Clear(){
121    //    //
122  };  };
123    
124  void ToFLevel2::Delete(){ //ELENA  void ToFLevel2::Delete(Option_t *t){ //ELENA
125    //    //
126    if(ToFTrk){    if(ToFTrk){
127        ToFTrk->Delete(); //ELENA        ToFTrk->Delete(); //ELENA
# Line 216  Int_t ToFLevel2::GetNHitPaddles(Int_t pl Line 221  Int_t ToFLevel2::GetNHitPaddles(Int_t pl
221      return npad;      return npad;
222  };  };
223    
224    //wm Nov 08
225  //gf Apr 07  //gf Apr 07
226  /**  /**
227   * 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:
228   * 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
229   *  and needs a revision.   * only the paddle hitted by the track gets a dEdx value and the other
230     * paddles are set to zero, the output is just the dEdx of the hitted
231     * paddle in each layer!
232     * The "adcfl" option is not very useful (an artificial dEdx is per
233     * definition= 1 mip and not a real measurement), anyway left in the code
234   * @param notrack Track Number   * @param notrack Track Number
235   * @param plane Plane index (0,1,2,3,4,5)   * @param plane Plane index (0,1,2,3,4,5)
236   * @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 390  void ToFLevel2::GetPMTIndex(Int_t ind, I
390    
391    
392    
393    //  wm Nov 08 revision - saturation values included
394  /// gf Apr 07  /// gf Apr 07
   
395  /**  /**
396   * Method to get the dEdx from a given ToF paddle.   * Method to get the dEdx from a given ToF paddle.
397     * If two PMTs are good, the mean dEdx of both PMTs is taken, otherwise
398     * just the dEdx of the "good" PMT. If both PMTs are above saturation => dEdx=1000
399   * @param notrack Track Number   * @param notrack Track Number
400   * @param Paddle index (0,1,...,23).   * @param Paddle index (0,1,...,23).
401   * @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 404  void ToFLevel2::GetPMTIndex(Int_t ind, I
404   */   */
405  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){
406    
407    /*
408    Float_t  PMTsat[48] = {
409    3162.14, 3165.48, 3153.85, 3085.73, 3089.65, 3107.64, 3097.52, 3078.37,
410    3130.05, 3087.07, 3112.22, 3102.92, 3080.58, 3092.55, 3087.94, 3125.03,
411    3094.09, 3143.16, 3125.51, 3181.27, 3092.09, 3124.98, 3069.3, 3095.53,
412    3097.11, 3133.53, 3114.73, 3113.01, 3091.19, 3097.99, 3033.84, 3134.98,
413    3081.37, 3111.04, 3066.77, 3108.17, 3133, 3111.06, 3052.52, 3140.66,
414    3106.33, 3094.85, 3150.85, 3118.8, 3096.24, 3118.47,3111.36, 3117.11 } ;
415    */
416    
417    // new values from Napoli dec 2008
418    Float_t  PMTsat[48] = {
419    3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
420    3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
421    3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
422    3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
423    3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
424    3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };
425    
426    for (Int_t i=0; i<48;i++) PMTsat[i] = PMTsat[i] - 5.;  // safety margin
427    
428    
429    PadEdx = 0.;    PadEdx = 0.;
430    SatWarning = 1000;  //  SatWarning = 1000;
431      SatWarning = 0;   // 0=good, increase for each bad PMT
432    
433    Float_t dEdx[48] = {0};    Float_t dEdx[48] = {0};
434    Int_t pmt_id = -1;    Int_t pmt_id = -1;
# Line 427  void ToFLevel2::GetdEdxPaddle(Int_t notr Line 460  void ToFLevel2::GetdEdxPaddle(Int_t notr
460        adcraw[pmtright] = pmt->adc;        adcraw[pmtright] = pmt->adc;
461      }      }
462    }    }
463    
464        
465    for (Int_t i=0; i<trk->npmtadc; i++){    for (Int_t i=0; i<trk->npmtadc; i++){
466    
# Line 439  void ToFLevel2::GetdEdxPaddle(Int_t notr Line 473  void ToFLevel2::GetdEdxPaddle(Int_t notr
473      }      }
474    }    }
475    
476    if( adcraw[pmtleft] >3000 || adcraw[pmtright] >3000)SatWarning=1;  
477        //  if( adcraw[pmtleft] >3000 || adcraw[pmtright] >3000)SatWarning=1;  //old version
478    if(dEdx[pmtleft]!=0 && dEdx[pmtright]!=0){  
479      PadEdx = (dEdx[pmtleft]+dEdx[pmtright])*0.5;  // Increase SatWarning Counter for each PMT>Sat
480    }    if( adcraw[pmtleft] > PMTsat[pmtleft])SatWarning++;  
481    if(dEdx[pmtleft]==0 && dEdx[pmtright]!=0){    if( adcraw[pmtright] > PMTsat[pmtright])SatWarning++;
482      PadEdx = dEdx[pmtright];  
483    }  // if ADC  > sat set dEdx=1000
484    if(dEdx[pmtleft]!=0 && dEdx[pmtright]==0){    if( adcraw[pmtleft] > PMTsat[pmtleft]) dEdx[pmtleft] = 1000.;
485      PadEdx = dEdx[pmtleft];    if( adcraw[pmtright] > PMTsat[pmtright]) dEdx[pmtright] = 1000. ;
486    }  
487    // if two PMT are good, take mean dEdx, otherwise only the good dEdx
488      if(dEdx[pmtleft]<1000 && dEdx[pmtright]<1000) PadEdx = (dEdx[pmtleft]+dEdx[pmtright])*0.5;
489      if(dEdx[pmtleft]==1000 && dEdx[pmtright]<1000) PadEdx = dEdx[pmtright];  
490      if(dEdx[pmtleft]<1000 && dEdx[pmtright]==1000) PadEdx = dEdx[pmtleft];
491        
   return;  
492  };  };
493  //  //
494    
# Line 508  TString ToFLevel2::GetPMTName(Int_t ind) Line 545  TString ToFLevel2::GetPMTName(Int_t ind)
545        
546  };  };
547    
548    // wm jun 08
 // gf Apr 07  
549  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){
550    return GetPaddleIdOfTrack(xtr ,ytr ,plane, 0.4);
551    }
552    
553    // gf Apr 07
554    Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane, Float_t margin){
555      
556    Double_t xt,yt,xl,xh,yl,yh;    Double_t xt,yt,xl,xh,yl,yh;
557        
558    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 582  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
582      yh =  33.0/2. ;      yh =  33.0/2. ;
583      if ((yt>yl)&&(yt<yh)) {      if ((yt>yl)&&(yt<yh)) {
584        for (Int_t i1=0; i1<8;i1++){        for (Int_t i1=0; i1<8;i1++){
585          xl = tof11_x[i1] - (5.1-0.4)/2. ;          xl = tof11_x[i1] - (5.1-margin)/2. ;
586          xh = tof11_x[i1] + (5.1-0.4)/2. ;          xh = tof11_x[i1] + (5.1-margin)/2. ;
587          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
588        }        }
589      }      }
# Line 559  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 600  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
600            
601      if ((xt>xl)&&(xt<xh)) {      if ((xt>xl)&&(xt<xh)) {
602        for (Int_t i1=0; i1<6;i1++){        for (Int_t i1=0; i1<6;i1++){
603          yl = tof12_y[i1] - (5.5-0.4)/2. ;          yl = tof12_y[i1] - (5.5-margin)/2. ;
604          yh = tof12_y[i1] + (5.5-0.4)/2. ;          yh = tof12_y[i1] + (5.5-margin)/2. ;
605          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;
606        }        }
607      }      }
# Line 577  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 618  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
618            
619      if ((xt>xl)&&(xt<xh)) {      if ((xt>xl)&&(xt<xh)) {
620        for (Int_t i1=0; i1<2;i1++){        for (Int_t i1=0; i1<2;i1++){
621          yl = tof21_y[i1] - (7.5-0.4)/2. ;          yl = tof21_y[i1] - (7.5-margin)/2. ;
622          yh = tof21_y[i1] + (7.5-0.4)/2. ;          yh = tof21_y[i1] + (7.5-margin)/2. ;
623          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;
624        }        }
625      }      }
# Line 594  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 635  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
635            
636      if ((yt>yl)&&(yt<yh)) {      if ((yt>yl)&&(yt<yh)) {
637        for (Int_t i1=0; i1<2;i1++){        for (Int_t i1=0; i1<2;i1++){
638          xl = tof22_x[i1] - (9.0-0.4)/2. ;          xl = tof22_x[i1] - (9.0-margin)/2. ;
639          xh = tof22_x[i1] + (9.0-0.4)/2. ;          xh = tof22_x[i1] + (9.0-margin)/2. ;
640          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
641        }        }
642      }      }
# Line 611  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 652  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
652            
653      if ((yt>yl)&&(yt<yh)) {      if ((yt>yl)&&(yt<yh)) {
654        for (Int_t i1=0; i1<3;i1++){        for (Int_t i1=0; i1<3;i1++){
655          xl = tof31_x[i1] - (6.0-0.4)/2. ;          xl = tof31_x[i1] - (6.0-margin)/2. ;
656          xh = tof31_x[i1] + (6.0-0.4)/2. ;          xh = tof31_x[i1] + (6.0-margin)/2. ;
657          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
658        }        }
659      }      }
# Line 628  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 669  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
669            
670      if ((xt>xl)&&(xt<xh)) {      if ((xt>xl)&&(xt<xh)) {
671        for (Int_t i1=0; i1<3;i1++){        for (Int_t i1=0; i1<3;i1++){
672          yl = tof32_y[i1] - (5.0-0.4)/2. ;          yl = tof32_y[i1] - (5.0-margin)/2. ;
673          yh = tof32_y[i1] + (5.0-0.4)/2. ;          yh = tof32_y[i1] + (5.0-margin)/2. ;
674          if ((yt>yl)&&(yt<yh)) paddleidoftrack=i1;          if ((yt>yl)&&(yt<yh)) paddleidoftrack=i1;
675        }        }
676      }      }
# Line 696  void ToFLevel2::GetPMTPaddle(Int_t pmt_i Line 737  void ToFLevel2::GetPMTPaddle(Int_t pmt_i
737  // gf Apr 07  // gf Apr 07
738    
739  void ToFLevel2::GetPaddlePMT(Int_t paddle, Int_t &pmtleft, Int_t &pmtright){  void ToFLevel2::GetPaddlePMT(Int_t paddle, Int_t &pmtleft, Int_t &pmtright){
740      pmtleft=paddle*2;
741    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;  
   }  
     
742    return;    return;
743  }  }
744    
# Line 930  void ToFLevel2::GetPaddleGeometry(Int_t Line 852  void ToFLevel2::GetPaddleGeometry(Int_t
852   */   */
853  Int_t ToFLevel2::GetPaddleid(Int_t plane, Int_t paddle)  Int_t ToFLevel2::GetPaddleid(Int_t plane, Int_t paddle)
854  {  {
   
855    Int_t padid=-1;    Int_t padid=-1;
856    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;  
   }  
857    
858    if(plane == 1){    int somma=0;
859      padid=pads11+paddle;    int np=plane;
860    }    for(Int_t j=0; j<np; j++){
861        somma+=pads[j];
   if(plane == 2){  
     padid=pads11+pads12+paddle;  
862    }    }
863      padid=paddle+somma;
   if(plane == 3){  
     padid=pads11+pads12+pads21+paddle;  
   }  
   
   if(plane == 4){  
     padid=pads11+pads12+pads21+pads22+paddle;  
   }  
   
   if(plane == 5){  
     padid=pads11+pads12+pads21+pads22+pads31+paddle;  
   }  
   
864    return padid;    return padid;
865    
866  }  }
# Line 994  void ToFLevel2::GetPaddlePlane(Int_t pad Line 891  void ToFLevel2::GetPaddlePlane(Int_t pad
891      return;      return;
892    }    }
893    
894    if(7<pad<14){    if((7<pad)&&(pad<14)){
895      plane=1;      plane=1;
896      paddle=pad-pads11;      paddle=pad-pads11;
897      return;      return;
898    }    }
899        
900    if(13<pad<16){    if((13<pad)&&(pad<16)){
901      plane=2;      plane=2;
902      paddle=pad-pads11-pads12;      paddle=pad-pads11-pads12;
903      return;      return;
904    }    }
905    
906    if(15<pad<18){    if((15<pad)&&(pad<18)){
907      plane=3;      plane=3;
908      paddle=pad-pads11-pads12-pads21;      paddle=pad-pads11-pads12-pads21;
909      return;      return;
910    }    }
911    
912    if(17<pad<21){    if((17<pad)&&(pad<21)){
913      plane=4;      plane=4;
914      paddle=pad-pads11-pads12-pads21-pads22;      paddle=pad-pads11-pads12-pads21-pads22;
915      return;      return;
916    }    }
917    
918    if(20<pad<24){    if((20<pad)&&(pad<24)){
919      plane=5;      plane=5;
920      paddle=pad-pads11-pads12-pads21-pads22-pads31;      paddle=pad-pads11-pads12-pads21-pads22-pads31;
921      return;      return;
# Line 1049  Int_t ToFLevel2::GetNPaddle(Int_t plane) Line 946  Int_t ToFLevel2::GetNPaddle(Int_t plane)
946    
947  }  }
948    
 ////////////////////////////////////////////////////  
949    
950    
951    /// wm feb 08
952    
953    /**
954     * Method to calculate Beta from the 12 single measurements
955     * we check the individual weights for artificial TDC values, then calculate
956     * am mean beta for the first time. In a second step we loop again through
957     * the single measurements, checking for the residual from the mean
958     * The cut on the residual reject measurements > "x"-sigma. A chi2 value is
959     * calculated, furthermore a "quality" value by adding the weights which
960     * are finally used. If all measurements are taken, "quality" will be = 22.47.
961     * A chi2 cut around 3-4 and a quality-cut > 20 is needed for clean beta
962     * measurements like antiprotons etc.
963     * The Level2 output is derived in the fortran routines using: 10.,10.,20.
964     * @param notrack Track Number
965     * @param cut on residual: difference between single measurement and mean
966     * @param cut on "quality"
967     * @param cut on chi2
968     */
969    
970    Float_t ToFLevel2::CalcBeta(Int_t notrack, Float_t resmax, Float_t qualitycut, Float_t chi2cut){
971    
972    //  cout<<" in CalcBeta "<<resmax<<" "<<chi2cut<<" "<<qualitycut<<endl;
973    
974      Float_t bxx = 100.;
975      //
976      ToFTrkVar *trk = GetToFTrkVar(notrack);
977      if(!trk) return 0; //ELENA
978    
979    
980      Float_t chi2,xhelp,beta_mean;
981      Float_t w_i[12],quality,sw,sxw,res,betachi,beta_mean_inv;
982      Float_t b[12],tdcfl;
983      Int_t  pmt_id,pmt_plane;
984    
985      for (Int_t i=0; i<12; i++){
986        b[i] = trk->beta[i];
987                                  }
988          
989    
990    //========================================================================
991    //---  Find out ToF layers with artificial TDC values & fill vector    ---
992    //========================================================================
993    
994    Float_t  w_il[6];
995    
996         for (Int_t jj=0; jj<6;jj++) {
997             w_il[jj] = 1000.;
998                                     }
999    
1000    
1001      for (Int_t i=0; i<trk->npmttdc; i++){
1002        //
1003        pmt_id = (trk->pmttdc).At(i);
1004        pmt_plane = GetPlaneIndex(pmt_id);
1005        tdcfl = (trk->tdcflag).At(i);
1006        if (w_il[pmt_plane] != 1.) w_il[pmt_plane] = tdcfl; //tdcflag
1007                                         };
1008      
1009    //========================================================================
1010    //---  Set weights for the 12 measurements using information for top and bottom:
1011    //---  if no measurements: weight = set to very high value=> not used
1012    //---  top or bottom artificial: weight*sqrt(2)
1013    //---  top and bottom artificial: weight*sqrt(2)*sqrt(2)
1014    //========================================================================
1015    
1016    Int_t itop[12] = {0,0,1,1,2,2,3,3,0,0,1,1};
1017    Int_t ibot[12] = {4,5,4,5,4,5,4,5,2,3,2,3};
1018    
1019         xhelp= 1E09;
1020      
1021         for (Int_t jj=0; jj<12;jj++) {
1022         if (jj<4)           xhelp = 0.11;    // S1-S3
1023         if ((jj>3)&&(jj<8)) xhelp = 0.18;    // S2-S3
1024         if (jj>7)           xhelp = 0.28;    // S1-S2
1025         if ((w_il[itop[jj]] == 1000.) && (w_il[ibot[jj]] == 1000.)) xhelp = 1E09;
1026         if ((w_il[itop[jj]] == 1) || (w_il[ibot[jj]] == 1.)) xhelp = xhelp*1.414 ;
1027         if ((w_il[itop[jj]] == 1) && (w_il[ibot[jj]] == 1.)) xhelp = xhelp*2. ;
1028    
1029         w_i[jj] = 1./xhelp;
1030                                      }
1031    
1032    
1033    //========================================================================
1034    //--- Calculate mean beta for the first time -----------------------------
1035    //--- We are using "1/beta" since its error is gaussian ------------------
1036    //========================================================================
1037    
1038          Int_t icount=0;
1039          sw=0.;
1040          sxw=0.;
1041          beta_mean=100.;
1042    
1043              for (Int_t jj=0; jj<12;jj++){
1044            if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.))
1045             {
1046                icount= icount+1;
1047                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1048                sw =sw + w_i[jj]*w_i[jj] ;
1049    
1050             }
1051             }
1052    
1053          if (icount>0) beta_mean=1./(sxw/sw);
1054          beta_mean_inv = 1./beta_mean;
1055    
1056    //========================================================================
1057    //--- Calculate beta for the second time, use residuals of the single
1058    //--- measurements to get a chi2 value
1059    //========================================================================
1060    
1061          icount=0;
1062          sw=0.;
1063          sxw=0.;
1064          betachi = 100.;
1065          chi2 = 0.;
1066          quality=0.;
1067    
1068    
1069              for (Int_t jj=0; jj<12;jj++){
1070           if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.)&&(w_i[jj]>0.01)) {
1071                res = beta_mean_inv - (1./b[jj]) ;
1072                if (fabs(res*w_i[jj])<resmax)          {;
1073                chi2 = chi2 + pow((res*w_i[jj]),2) ;
1074                icount= icount+1;
1075                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1076                sw =sw + w_i[jj]*w_i[jj] ;
1077                                                   }
1078                                                                            }
1079                                          }
1080          quality = sqrt(sw) ;
1081    
1082          if (icount==0) chi2 = 1000.;
1083          if (icount>0) chi2 = chi2/(icount) ;
1084          if (icount>0) betachi=1./(sxw/sw);
1085    
1086       bxx = 100.;
1087       if ((chi2 < chi2cut)&&(quality>qualitycut)) bxx = betachi;
1088      //
1089      return(bxx);
1090    };
1091    
1092    
1093    ////////////////////////////////////////////////////
1094    ////////////////////////////////////////////////////
1095    
1096    
1097  /**  /**
1098   * 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 1140  void ToFLevel2::GetLevel2Struct(cToFLeve
1140        }        }
1141    } //ELENA    } //ELENA
1142  }  }
1143    
1144    
1145    //
1146    // Reprocessing tool // Emiliano 08/04/07
1147    //
1148    Int_t ToFLevel2::Process(TrkLevel2 *trk, TrigLevel2 *trg, GL_RUN *run, OrbitalInfo *orb, Bool_t force){
1149      //
1150      // Copiare qui qualcosa di simile a calonuclei per evitare di riprocessare sempre tutto
1151      //
1152    
1153    
1154    
1155    
1156      //
1157      // structures to communicate with F77
1158      //
1159      extern struct ToFInput  tofinput_;
1160      extern struct ToFOutput tofoutput_;
1161      //
1162      // DB connection
1163      //
1164      TString host;
1165      TString user;
1166      TString psw;
1167      const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
1168      const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
1169      const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
1170      if ( !pamdbhost ) pamdbhost = "";
1171      if ( !pamdbuser ) pamdbuser = "";
1172      if ( !pamdbpsw ) pamdbpsw = "";
1173      if ( strcmp(pamdbhost,"") ) host = pamdbhost;
1174      if ( strcmp(pamdbuser,"") ) user = pamdbuser;
1175      if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
1176      //
1177      //
1178      TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
1179      if ( !dbc->IsConnected() ) return 1;
1180      stringstream myquery;
1181      myquery.str("");
1182      myquery << "SET time_zone='+0:00'";
1183      dbc->Query(myquery.str().c_str());
1184      GL_PARAM *glparam = new GL_PARAM();
1185      glparam->Query_GL_PARAM(1,1,dbc); // parameters stored in DB in GL_PRAM table
1186      trk->LoadField(glparam->PATH+glparam->NAME);
1187      //
1188      Bool_t defcal = true;
1189      Int_t error=glparam->Query_GL_PARAM(run->RUNHEADER_TIME,201,dbc); // parameters stored in DB in GL_PRAM table
1190      if ( error<0 ) {
1191        return(1);
1192      };
1193      printf(" Reading ToF parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
1194      if ( (UInt_t)glparam->TO_TIME != (UInt_t)4294967295UL ) defcal = false;
1195      //
1196      Int_t nlen = (Int_t)(glparam->PATH+glparam->NAME).Length();
1197      rdtofcal((char *)(glparam->PATH+glparam->NAME).Data(),&nlen);
1198      //
1199      Int_t adc[4][12];
1200      Int_t tdc[4][12];
1201      Float_t tdcc[4][12];
1202      //
1203      // process tof data
1204      //
1205      for (Int_t hh=0; hh<12;hh++){
1206        for (Int_t kk=0; kk<4;kk++){
1207               adc[kk][hh] = 4095;
1208               tdc[kk][hh] = 4095;
1209               tdcc[kk][hh] = 4095.;
1210               tofinput_.adc[hh][kk] = 4095;
1211               tofinput_.tdc[hh][kk] = 4095;
1212        };
1213      };
1214      Int_t ntrkentry = 0;
1215      Int_t npmtentry = 0;
1216      Int_t gg = 0;
1217      Int_t hh = 0;
1218      Int_t adcf[48];
1219      memset(adcf, 0, 48*sizeof(Int_t));
1220      Int_t tdcf[48];
1221      memset(tdcf, 0, 48*sizeof(Int_t));
1222      for (Int_t pm=0; pm < this->ntrk() ; pm++){
1223         ToFTrkVar *ttf = this->GetToFTrkVar(pm);
1224         for ( Int_t nc=0; nc < ttf->npmttdc; nc++){
1225                if ( (ttf->tdcflag).At(nc) != 0 ) tdcf[(ttf->pmttdc).At(nc)] = 1;
1226         };
1227         for ( Int_t nc=0; nc < ttf->npmtadc; nc++){
1228                if ( (ttf->adcflag).At(nc) != 0 ) adcf[(ttf->pmtadc).At(nc)] = 1;
1229         };
1230      };
1231      //
1232      for (Int_t pm=0; pm < this->npmt() ; pm++){
1233         ToFPMT *pmt = this->GetToFPMT(pm);
1234         this->GetPMTIndex(pmt->pmt_id, gg, hh);
1235         if ( adcf[pmt->pmt_id] == 0 ){
1236                 tofinput_.adc[gg][hh] = (int)pmt->adc;
1237                 adc[hh][gg] = (int)pmt->adc;
1238         };
1239         if ( tdcf[pmt->pmt_id] == 0 ){
1240                 tofinput_.tdc[gg][hh] = (int)pmt->tdc;
1241                 tdc[hh][gg] = (int)pmt->tdc;
1242         };
1243         tdcc[hh][gg] = (float)pmt->tdc_tw;
1244         // Int_t pppid = this->GetPMTid(hh,gg);
1245         //      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);
1246      };
1247      //
1248      Int_t unpackError = this->unpackError;
1249      //
1250      for (Int_t hh=0; hh<5;hh++){
1251         tofinput_.patterntrig[hh]=trg->patterntrig[hh];
1252      };
1253      //
1254      this->Clear();
1255      //
1256          Int_t pmt_id = 0;
1257          ToFPMT *t_pmt = new ToFPMT();
1258          if(!(this->PMT)) this->PMT = new TClonesArray("ToFPMT",12); //ELENA
1259          TClonesArray &tpmt = *this->PMT;
1260          ToFTrkVar *t_tof = new ToFTrkVar();
1261          if(!(this->ToFTrk)) this->ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
1262          TClonesArray &t = *this->ToFTrk;
1263          //
1264          //
1265          // Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract common and track-related  variables.
1266          //
1267          npmtentry = 0;
1268          //
1269          ntrkentry = 0;
1270          //
1271          // Calculate tracks informations from ToF alone
1272          //
1273          tofl2com();
1274          //
1275          memcpy(this->tof_j_flag,tofoutput_.tof_j_flag,6*sizeof(Int_t));
1276          //
1277          t_tof->trkseqno = -1;
1278          //
1279          // and now we must copy from the output structure to the level2 class:
1280          //
1281          t_tof->npmttdc = 0;
1282          //
1283          for (Int_t hh=0; hh<12;hh++){
1284            for (Int_t kk=0; kk<4;kk++){
1285              if ( tofoutput_.tofmask[hh][kk] != 0 ){
1286                pmt_id = this->GetPMTid(kk,hh);
1287                t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
1288                t_tof->tdcflag.AddAt(tofoutput_.tdcflagtof[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
1289                t_tof->npmttdc++;
1290              };
1291            };
1292          };
1293          for (Int_t kk=0; kk<13;kk++){
1294            t_tof->beta[kk] = tofoutput_.betatof_a[kk];
1295          }
1296          //
1297          t_tof->npmtadc = 0;
1298          for (Int_t hh=0; hh<12;hh++){
1299            for (Int_t kk=0; kk<4;kk++){
1300              if ( tofoutput_.adctof_c[hh][kk] < 1000 ){
1301                t_tof->dedx.AddAt(tofoutput_.adctof_c[hh][kk],t_tof->npmtadc);
1302                pmt_id = this->GetPMTid(kk,hh);
1303                t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
1304                t_tof->adcflag.AddAt(tofoutput_.adcflagtof[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
1305                t_tof->npmtadc++;
1306              };
1307            };
1308          };
1309          //
1310          memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
1311          memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
1312          memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
1313          memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
1314          //
1315          new(t[ntrkentry]) ToFTrkVar(*t_tof);
1316          ntrkentry++;
1317          t_tof->Clear();
1318          //
1319          //
1320          //
1321          t_pmt->Clear();
1322          //
1323          for (Int_t hh=0; hh<12;hh++){
1324            for (Int_t kk=0; kk<4;kk++){
1325             // new WM
1326              if ( tofoutput_.tdc_c[hh][kk] < 4095 || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){
1327    //          if ( tdcc[kk][hh] < 4095. || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){
1328                //
1329                t_pmt->pmt_id = this->GetPMTid(kk,hh);
1330                t_pmt->tdc_tw = tofoutput_.tdc_c[hh][kk];
1331                t_pmt->adc = (Float_t)adc[kk][hh];
1332                t_pmt->tdc = (Float_t)tdc[kk][hh];
1333                //
1334                new(tpmt[npmtentry]) ToFPMT(*t_pmt);
1335                npmtentry++;
1336                t_pmt->Clear();
1337              };
1338            };
1339          };
1340          //
1341          // Calculate track-related variables
1342          //
1343          if ( trk->ntrk() > 0 ){
1344            //
1345            // We have at least one track
1346            //
1347            //
1348            // Run over tracks
1349            //
1350            for(Int_t nt=0; nt < trk->ntrk(); nt++){
1351              //
1352              TrkTrack *ptt = trk->GetStoredTrack(nt);
1353              //
1354              // Copy the alpha vector in the input structure
1355              //
1356              for (Int_t e = 0; e < 5 ; e++){
1357                tofinput_.al_pp[e] = ptt->al[e];
1358              };
1359              //
1360              // Get tracker related variables for this track
1361              //
1362              toftrk();
1363              //
1364              // Copy values in the class from the structure (we need to use a temporary class to store variables).
1365              //
1366              t_tof->npmttdc = 0;
1367              for (Int_t hh=0; hh<12;hh++){
1368                for (Int_t kk=0; kk<4;kk++){
1369                  if ( tofoutput_.tofmask[hh][kk] != 0 ){
1370                    pmt_id = this->GetPMTid(kk,hh);
1371                    t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
1372                    t_tof->tdcflag.AddAt(tofoutput_.tdcflag[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
1373                    t_tof->npmttdc++;
1374                  };
1375                };
1376              };
1377              for (Int_t kk=0; kk<13;kk++){
1378                t_tof->beta[kk] = tofoutput_.beta_a[kk];
1379              };
1380              //
1381              t_tof->npmtadc = 0;
1382              for (Int_t hh=0; hh<12;hh++){
1383                for (Int_t kk=0; kk<4;kk++){
1384                  if ( tofoutput_.adc_c[hh][kk] < 1000 ){
1385                    t_tof->dedx.AddAt(tofoutput_.adc_c[hh][kk],t_tof->npmtadc);
1386                    pmt_id = this->GetPMTid(kk,hh);
1387                    t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
1388                    t_tof->adcflag.AddAt(tofoutput_.adcflag[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
1389                    t_tof->npmtadc++;
1390                  };
1391                };
1392              };
1393              //
1394              memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
1395              memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
1396              memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
1397              memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
1398              //
1399              // Store the tracker track number in order to be sure to have shyncronized data during analysis
1400              //
1401              t_tof->trkseqno = nt;
1402              //
1403              // create a new object for this event with track-related variables
1404              //
1405              new(t[ntrkentry]) ToFTrkVar(*t_tof);
1406              ntrkentry++;
1407              t_tof->Clear();
1408              //
1409            }; // loop on all the tracks
1410          //
1411          this->unpackError = unpackError;
1412          if ( defcal ){
1413            this->default_calib = 1;
1414          } else {
1415            this->default_calib = 0;
1416          };
1417     };
1418    
1419    
1420    
1421      return(0);
1422    }
1423    
1424    
1425    
1426    
1427    
1428    
1429    
1430    
1431    
1432    
1433    
1434    
1435    
1436    
1437    
1438    
1439    
1440    
1441    
1442    
1443    
1444    
1445    
1446    
1447    ToFdEdx::ToFdEdx()
1448    {
1449      memset(conn,0,12*sizeof(Bool_t));
1450      memset(ts,0,12*sizeof(UInt_t));
1451      memset(te,0,12*sizeof(UInt_t));
1452      Define_PMTsat();
1453      Clear();
1454    }
1455    //------------------------------------------------------------------------
1456    void ToFdEdx::CheckConnectors(UInt_t atime, GL_PARAM *glparam, TSQLServer *dbc)
1457    {
1458      for(int i=0; i<12; i++){
1459        if(atime<=ts[i] || atime>te[i]){
1460          Int_t error=glparam->Query_GL_PARAM(atime,210+i,dbc); // parameters stored in DB in GL_PRAM table
1461          if ( error<0 ) {
1462            conn[i]=false;
1463            ts[i]=0;
1464            te[i]=numeric_limits<UInt_t>::max();
1465          };
1466          if ( !error ){
1467            conn[i]=true;
1468            ts[i]=glparam->FROM_TIME;
1469            te[i]=glparam->TO_TIME;
1470          }
1471          if ( error>0 ){
1472            conn[i]=false;
1473            ts[i]=glparam->TO_TIME;
1474            TSQLResult *pResult;
1475            TSQLRow *row;
1476            TString query= Form("SELECT FROM_TIME FROM GL_PARAM WHERE TYPE=%i AND FROM_TIME>=%i ORDER BY FROM_TIME ASC LIMIT 1;",210+i,atime);
1477            pResult=dbc->Query(query.Data());
1478            if(!pResult->GetRowCount()){
1479              te[i]=numeric_limits<UInt_t>::max();
1480            }else{
1481              row=pResult->Next();
1482              te[i]=(UInt_t)atoll(row->GetField(0));
1483            }
1484          }
1485          //
1486          
1487        }
1488      }
1489    
1490    }
1491    //------------------------------------------------------------------------
1492    void ToFdEdx::Clear(Option_t *option)
1493    {
1494      //
1495      // Set arrays and initialize structure
1496    
1497      eDEDXpmt.Set(48);    eDEDXpmt.Reset(-1);   // Set array size  and reset structure
1498      eZpmt.Set(48);       eZpmt.Reset(-1);
1499      eDEDXpad.Set(24);    eDEDXpad.Reset(-1);
1500      eZpad.Set(24);       eZpad.Reset(-1);
1501      eDEDXlayer.Set(6);   eDEDXlayer.Reset(-1);
1502      eZlayer.Set(6);      eZlayer.Reset(-1);
1503      eDEDXplane.Set(3);   eDEDXplane.Reset(-1);
1504      eZplane.Set(3);      eZplane.Reset(-1);
1505      INFOpmt.Set(48);     INFOpmt.Reset(0);
1506      INFOlayer.Set(6);    INFOlayer.Reset(0);
1507      //
1508    };
1509    
1510    //------------------------------------------------------------------------
1511    void ToFdEdx::Print(Option_t *option)
1512    {
1513      //
1514      printf("========================================================================\n");
1515    
1516    };
1517    
1518    
1519    //------------------------------------------------------------------------
1520    // void ToFdEdx::InitPar(TString parname, TString parfile)
1521    // {
1522    //   // expensive function - call it once/run
1523    
1524    
1525    
1526    //   ReadParAtt(            Form("%s/attenuation.txt"              , pardir) );
1527    //   ReadParPos(            Form("%s/desaturation_position.txt"    , pardir) );
1528    //   ReadParBBneg(          Form("%s/BetheBloch.txt"               , pardir) );
1529    //   ReadParBBpos(          Form("%s/BetheBloch_betagt1.txt"       , pardir) );
1530    //   ReadParDesatBB(        Form("%s/desaturation_beta.txt"        , pardir) );
1531    
1532    // };
1533    
1534    
1535    //------------------------------------------------------------------------
1536    void ToFdEdx::Process(UInt_t atime, Float_t betamean, Float_t *xtr_tof, Float_t *ytr_tof, pamela::tof::TofEvent *tofl0 )
1537    {
1538      // the parameters should be already initialised by InitPar()
1539    
1540    
1541      Clear();
1542    
1543    
1544    
1545      //  Float_t betamean = fabs(trackTRK->GetToFTrack()->beta[12]);
1546    
1547      if(betamean<0.05 || betamean>2){
1548        for(int i=0;i<48;i++)INFOpmt[i]=1;
1549      }
1550    
1551     // define angle:  
1552      double dx   = xtr_tof[1] - xtr_tof[5];
1553      double dy   = ytr_tof[0] - ytr_tof[4];
1554      double dr   = sqrt(dx*dx+dy*dy);
1555      double theta=atan(dr/76.81);
1556    
1557    
1558    
1559      //  TArrayF adc;
1560      Float_t adc[48];
1561    
1562      ToFLevel2 tf;
1563    
1564      for (Int_t gg=0; gg<4;gg++){
1565        for (Int_t hh=0; hh<12;hh++){
1566          //          tofinput_.tdc[hh][gg]=tofEvent->tdc[gg][hh];          
1567          int mm = tf.GetPMTid(gg,hh);        
1568          adc[mm]=tofl0->adc[gg][hh];
1569        };      
1570      };
1571      
1572      
1573      
1574      
1575      for( int ii=0; ii<48; ii++ ) {
1576        if( adc[ii] >= PMTsat[ii]-5 )  continue;
1577        if( adc[ii] <= 0. )            continue;
1578      
1579        double adcpC   = f_adcPC( adc[ii] );    // - adc conversion in pC
1580        double adccorr = adcpC*fabs(cos(theta));
1581    
1582             if(adccorr<=0.)           continue;
1583    
1584    
1585    
1586        //--------------------- TABLE OF PERIODS WITH HV PROBLEMS ----------------------------
1587    
1588        int Aconn=conn[0];    // PMT 0,20,22,24
1589        int Bconn=conn[1];    // PMT 6,12,26,34
1590        int Cconn=conn[2];    // PMT 4,14,28,32
1591        int Dconn=conn[3];    // PMT 2,8,10,30
1592        int Econn=conn[4];    // PMT 42,43,44,47
1593        int Fconn=conn[5];    // PMT 7,19,23,27
1594        int Gconn=conn[6];    // PMT 3,11,25,33
1595        int Hconn=conn[7];    // PMT 1,9,13,21
1596        int Iconn=conn[8];    // PMT 5,29,31,35
1597        int Lconn=conn[9];    // PMT 37,40,45,46
1598        int Mconn=conn[10];    // PMT 15,16,17,18
1599        int Nconn=conn[11];    // PMT 36,38,39,41
1600    
1601        //    int standard=0;
1602        if( false ) cout << Gconn << Iconn << Lconn <<endl;
1603        int S115B_ok=0;
1604        int S115B_break=0;
1605    
1606    //   if(atime>=1153660001 && atime<=1154375000)Dconn=1;
1607    //     else if(atime>=1155850001 && atime<=1156280000){
1608    //       Hconn=1;
1609    //       Nconn=1;
1610    //     }
1611    
1612    //  else if(atime>=1168490001 && atime<=1168940000)Dconn=1;
1613    //     else if(atime>=1168940001 && atime<=1169580000){
1614    //       Fconn=1;
1615    //       Mconn=1;
1616    //     }
1617    
1618    //  else if(atime>=1174665001 && atime<=1175000000)Bconn=1;
1619    //     else if(atime>=1176120001 && atime<=1176800000)Hconn=1;
1620    //     else if(atime>=1176800001 && atime<=1178330000)Econn=1;
1621    //     else if(atime>=1178330001 && atime<=1181322000)Hconn=1;
1622    //     else if(atime>=1182100001 && atime<=1183030000)Aconn=1;
1623    //     else if(atime>=1184000001 && atime<=1184570000)Hconn=1;
1624    //     else if(atime>=1185090001 && atime<=1185212000)Dconn=1;
1625    //     else if(atime>=1191100001 && atime<=1191940000)Dconn=1;
1626    //     else if(atime>=1196230001 && atime<=1196280000)Hconn=1;
1627    //     else if(atime>=1206100001 && atime<=1206375600)Cconn=1;
1628    //     else if(atime>=1217989201 && atime<=1218547800)Econn=1;
1629    //     else if(atime>=1225789201 && atime<=1226566800)Econn=1;
1630    //     else if(atime>=1229400901 && atime<=1229700000)Econn=1;
1631    //     else if(atime>=1230318001 && atime<=1230415200)Econn=1;
1632    //     else {
1633    //       standard=1;
1634    //     }
1635        if(atime<1158720000)S115B_ok=1;
1636        else S115B_break=1;
1637    
1638    
1639     //------------------------------------------------------------------------
1640    
1641    //---------------------------------------------------- Z reconstruction
1642    
1643    double adcHe, adcnorm, adclin, dEdx, Zeta;
1644    
1645     adcHe=-2;
1646     adcnorm=-2;
1647     adclin=-2;
1648     dEdx=-2;
1649     Zeta=-2;
1650    
1651    
1652    //  float ZetaH=-2;
1653    //  float dEdxH=-2;
1654    
1655    //  double day = (atime-1150000000)/84600;
1656    
1657        if(Aconn==1 && (ii==0 || ii==20 || ii==22 || ii==24)){
1658           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.675;
1659        }
1660        else if(Bconn==1 && (ii==6 || ii==12 || ii==26 || ii==34)){
1661           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/2.482;
1662        }
1663        else if(Cconn==1 && (ii==4 || ii==14 || ii==28 || ii==32)){
1664          adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.464;
1665        }
1666        else if(Dconn==1 && (ii==2 || ii==8 || ii==10 || ii==30)){
1667           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.995;
1668        }
1669        else if(Econn==1 && (ii==42 || ii==43 || ii==44 || ii==47)){
1670           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.273;
1671        }
1672        else if(Fconn==1 && (ii==7 || ii==19 || ii==23 || ii==27)){
1673           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.565;
1674        }
1675        else if(Mconn==1 && (ii==15 || ii==16 || ii==17 || ii==18)){
1676           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.565;
1677        }
1678        else if(Nconn==1 && (ii==36 || ii==38 || ii==39 || ii==41)){
1679           adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.018;
1680        }
1681        else if(Hconn==1 && (ii==1 || ii==13 || ii==21 || (ii==9&&S115B_ok==1))){
1682          adcHe   = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.84;
1683        }
1684        else if(S115B_break==1 && ii==9 && Hconn==0){
1685           adcHe   = f_att5B( ytr_tof[0] );   //N.B.: this function refers to the Carbon!!!
1686        }
1687        else if(S115B_break==1 && ii==9 && Hconn==1){
1688           adcHe   = (f_att5B( ytr_tof[0] ))/1.64;
1689        }
1690        else  adcHe   = Get_adc_he(ii, xtr_tof, ytr_tof);
1691    
1692        if(adcHe<=0)   continue;
1693    
1694        if(ii==9 && S115B_break==1)  adcnorm = f_pos5B(adccorr);
1695        else adcnorm = f_pos( (parPos[ii]), adccorr);
1696    
1697        if(adcnorm<=0) continue;
1698    
1699        if(ii==9 && S115B_break==1)  adclin  = 36.*adcnorm/adcHe;
1700        else  adclin  = 4.*adcnorm/adcHe;
1701    
1702        if(adclin<=0)  continue;
1703    
1704        double dEdxHe=-2;
1705        if(ii==9 && S115B_break==1){
1706          if( betamean <1. ) dEdxHe = f_BB5B( betamean );
1707          else                       dEdxHe = 33;
1708        } else {
1709          if( betamean <1. ) dEdxHe = f_BB( (parBBneg[ii]), betamean );
1710          else                       dEdxHe = parBBpos[ii];
1711        }
1712    
1713        if(dEdxHe<=0)  continue;
1714    
1715        if(ii==9 && S115B_break==1)  dEdx = f_desatBB5B( adclin );
1716        else  dEdx = f_desatBB((parDesatBB[ii]), adclin );
1717    
1718        if(dEdx<=0)    continue;
1719    
1720        if(ii==9 && S115B_break==1)  Zeta = sqrt(36.*(dEdx/dEdxHe));
1721        else  Zeta = sqrt(4.*(dEdx/dEdxHe));
1722    
1723        if(Zeta<=0)    continue;
1724    
1725    //--------------------- TIME DEPENDENCE ----------------------------------
1726    
1727    
1728    //      TArrayF &binx = TDx[ii];
1729    //      TArrayF &biny = TDy[ii];
1730    //      for(int k=0; k<200; k++) {
1731    //        if (day > binx[k]-2.5 && day<=binx[k]+2.5 && biny[k]>0)  {
1732    //       ZetaH=Zeta/biny[k]*6;
1733    //       dEdxH=dEdx/(pow(biny[k],2))*36;
1734    //        }
1735    //      }
1736    
1737    //      if(ZetaH!=-2)eZpmt[ii]=(Float_t)ZetaH;
1738    //      else eZpmt[ii]=(Float_t)Zeta;
1739    
1740    //      if(dEdxH!=-2)eDEDXpmt[ii]=(Float_t)dEdxH;
1741    //      else eDEDXpmt[ii]=(Float_t)dEdx;
1742    
1743    //      printf("%5d %8.2f %8.2f %8.2f  %8.2f %8.2f  %8.2f %5.4f \n",               ii, adcpC,  adccorr, adcHe, dEdxHe, dEdx, Zeta, betamean );
1744    
1745        eZpmt[ii]=(Float_t)Zeta;
1746        eDEDXpmt[ii]=(Float_t)dEdx;
1747    
1748    
1749     }  //end loop on 48 PMT
1750    
1751    //---------------------------------------------------  paddle + layer --------------------
1752    
1753      for(int j=0;j<48;j++){
1754        int k=100;
1755        if(j%2==0 || j==0)k=j/2;
1756        
1757        double zpdl=-1;
1758        
1759        if((j%2==0 || j==0) && eZpmt[j]!=-1 && eZpmt[j+1]!=-1){
1760          zpdl=0.5*(eZpmt[j]+eZpmt[j+1]);
1761        }else if((j%2==0 || j==0) && eZpmt[j]!=-1 && eZpmt[j+1]==-1){
1762          zpdl=eZpmt[j];
1763        }else if((j%2==0 || j==0) && eZpmt[j]==-1 && eZpmt[j+1]!=-1){
1764          zpdl=eZpmt[j+1];
1765        }
1766        
1767        if(j%2==0 || j==0)eZpad[k]= (Float_t)zpdl;
1768        
1769        if((j%2==0 || j==0)&&eZpad[k]!=-1){
1770          if(k>=0&&k<8)eZlayer[0]=eZpad[k];
1771          if(k>=8&&k<14)eZlayer[1]=eZpad[k];
1772          if(k>=14&&k<16)eZlayer[2]=eZpad[k];
1773          if(k>=16&&k<18)eZlayer[3]=eZpad[k];
1774          if(k>=18&&k<21)eZlayer[4]=eZpad[k];
1775          if(k>=21)eZlayer[5]=eZpad[k];
1776        }
1777    
1778        if(eZlayer[0]!=-1&&eZlayer[1]!=-1&&fabs(eZlayer[0]-eZlayer[1])<1.5)eZplane[0]=0.5*(eZlayer[0]+eZlayer[1]);
1779        else if(eZlayer[0]!=-1&&eZlayer[1]==-1)eZplane[0]=eZlayer[0];
1780        else if(eZlayer[1]!=-1&&eZlayer[0]==-1)eZplane[0]=eZlayer[1];
1781    
1782        if(eZlayer[2]!=-1&&eZlayer[3]!=-1&&fabs(eZlayer[2]-eZlayer[3])<1.5)eZplane[1]=0.5*(eZlayer[2]+eZlayer[3]);
1783        else if(eZlayer[2]!=-1&&eZlayer[3]==-1)eZplane[1]=eZlayer[2];
1784        else if(eZlayer[3]!=-1&&eZlayer[2]==-1)eZplane[1]=eZlayer[3];
1785    
1786        if(eZlayer[4]!=-1&&eZlayer[5]!=-1&&fabs(eZlayer[4]-eZlayer[5])<1.5)eZplane[2]=0.5*(eZlayer[4]+eZlayer[5]);
1787        else if(eZlayer[4]!=-1&&eZlayer[5]==-1)eZplane[2]=eZlayer[4];
1788        else if(eZlayer[5]!=-1&&eZlayer[4]==-1)eZplane[2]=eZlayer[5];
1789    
1790      }
1791    
1792      for(int jj=0;jj<48;jj++){
1793        int k=100;
1794        if(jj%2==0 || jj==0)k=jj/2;
1795        
1796        double dedxpdl=-1;
1797        
1798        if((jj%2==0 || jj==0) && eDEDXpmt[jj]!=-1 && eDEDXpmt[jj+1]!=-1){
1799          dedxpdl=0.5*(eDEDXpmt[jj]+eDEDXpmt[jj+1]);
1800        }else if((jj%2==0 || jj==0) && eDEDXpmt[jj]!=-1 && eDEDXpmt[jj+1]==-1){
1801          dedxpdl=eDEDXpmt[jj];
1802        }else if((jj%2==0 || jj==0) && eDEDXpmt[jj]==-1 && eDEDXpmt[jj+1]!=-1){
1803          dedxpdl=eDEDXpmt[jj+1];
1804        }
1805        
1806        if(jj%2==0 || jj==0)eDEDXpad[k]= (Float_t)dedxpdl;
1807        
1808        if((jj%2==0 || jj==0)&&eDEDXpad[k]!=-1){
1809          if(k>=0&&k<8)eDEDXlayer[0]=eDEDXpad[k];
1810          if(k>=8&&k<14)eDEDXlayer[1]=eDEDXpad[k];
1811          if(k>=14&&k<16)eDEDXlayer[2]=eDEDXpad[k];
1812          if(k>=16&&k<18)eDEDXlayer[3]=eDEDXpad[k];
1813          if(k>=18&&k<21)eDEDXlayer[4]=eDEDXpad[k];
1814          if(k>=21)eDEDXlayer[5]=eDEDXpad[k];
1815        }
1816    
1817        if(eDEDXlayer[0]!=-1&&eDEDXlayer[1]!=-1&&fabs(eDEDXlayer[0]-eDEDXlayer[1])<10)eDEDXplane[0]=0.5*(eDEDXlayer[0]+eDEDXlayer[1]);
1818        else if(eDEDXlayer[0]!=-1&&eDEDXlayer[1]==-1)eDEDXplane[0]=eDEDXlayer[0];
1819        else if(eDEDXlayer[1]!=-1&&eDEDXlayer[0]==-1)eDEDXplane[0]=eDEDXlayer[1];
1820    
1821        if(eDEDXlayer[2]!=-1&&eDEDXlayer[3]!=-1&&fabs(eDEDXlayer[2]-eDEDXlayer[3])<10)eDEDXplane[1]=0.5*(eDEDXlayer[2]+eDEDXlayer[3]);
1822        else if(eDEDXlayer[2]!=-1&&eDEDXlayer[3]==-1)eDEDXplane[1]=eDEDXlayer[2];
1823        else if(eDEDXlayer[3]!=-1&&eDEDXlayer[2]==-1)eDEDXplane[1]=eDEDXlayer[3];
1824    
1825        if(eDEDXlayer[4]!=-1&&eDEDXlayer[5]!=-1&&fabs(eDEDXlayer[4]-eDEDXlayer[5])<10)eDEDXplane[2]=0.5*(eDEDXlayer[4]+eDEDXlayer[5]);
1826        else if(eDEDXlayer[4]!=-1&&eDEDXlayer[5]==-1)eDEDXplane[2]=eDEDXlayer[4];
1827        else if(eDEDXlayer[5]!=-1&&eDEDXlayer[4]==-1)eDEDXplane[2]=eDEDXlayer[5];
1828    
1829      }
1830      
1831    
1832    
1833    };
1834    
1835    
1836    //------------------------------------------------------------------------
1837    void ToFdEdx::PrintTD()
1838    {
1839      for(int i=0; i<48; i++) {  
1840        TArrayF &binx = TDx[i];
1841        TArrayF &biny = TDy[i];
1842        for(int k=0; k<200; k++) {  // bin temporali
1843          printf("%d %d %f %f", i,k, binx[k], biny[k]);
1844          
1845        }
1846      }
1847    }
1848    
1849    
1850    //------------------------------------------------------------------------
1851    void ToFdEdx::Define_PMTsat()
1852    {
1853      Float_t  sat[48] = {
1854        3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
1855        3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
1856        3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
1857        3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
1858        3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
1859        3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };
1860      PMTsat.Set(48,sat);
1861    }
1862    
1863    //------------------------------------------------------------------------
1864    // void ToFdEdx::ReadParTD( Int_t ipmt, const char *fname )
1865    // {
1866    //   printf("read %s\n",fname);
1867    //   if(ipmt<0)  return;
1868    //   if(ipmt>47) return;
1869    //   FILE *fattin = fopen( fname , "r" );
1870    //   Float_t yTD[200],xTD[200];
1871    //   for(int j=0;j<200;j++){
1872    //     float x,y,ym,e;
1873    //     if(fscanf(fattin,"%f %f %f %f",
1874    //            &x, &y, &ym, &e )!=4) break;
1875    //     xTD[j]=x;
1876    //     if(ym>0&&fabs(y-ym)>1)  yTD[j]=ym;
1877    //     else                    yTD[j]=y;
1878    //   }
1879    //   TDx[ipmt].Set(200,xTD);
1880    //   TDy[ipmt].Set(200,yTD);
1881    //   fclose(fattin);
1882    // }
1883    
1884    //------------------------------------------------------------------------
1885    void ToFdEdx::ReadParBBpos( const char *fname )
1886    {
1887      printf("read %s\n",fname);
1888      parBBpos.Set(48);
1889      FILE *fattin = fopen( fname , "r" );
1890      for (int i=0; i<48; i++) {
1891        int   tid=0;
1892        float  tp;
1893        if(fscanf(fattin,"%d %f",
1894                  &tid, &tp )!=2) break;
1895        parBBpos[i]=tp;
1896      }
1897      fclose(fattin);
1898    }
1899    
1900    //------------------------------------------------------------------------
1901    void ToFdEdx::ReadParDesatBB( const char *fname )
1902    {
1903      printf("read %s\n",fname);
1904      FILE *fattin = fopen( fname , "r" );
1905      for (int i=0; i<48; i++) {
1906        int   tid=0;
1907        float  tp[3];
1908        if(fscanf(fattin,"%d %f %f %f",
1909                  &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
1910        parDesatBB[i].Set(3,tp);
1911      }
1912      fclose(fattin);
1913    }
1914    
1915    
1916    //------------------------------------------------------------------------
1917    void ToFdEdx::ReadParBBneg( const char *fname )
1918    
1919    {
1920      printf("read %s\n",fname);
1921      FILE *fattin = fopen( fname , "r" );
1922      for (int i=0; i<48; i++) {
1923        int   tid=0;
1924        float  tp[3];
1925        if(fscanf(fattin,"%d %f %f %f",
1926                  &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
1927        parBBneg[i].Set(3,tp);
1928      }
1929      fclose(fattin);
1930    }
1931    
1932    //------------------------------------------------------------------------
1933    void ToFdEdx::ReadParPos( const char *fname )
1934    {
1935      printf("read %s\n",fname);
1936      FILE *fattin = fopen( fname , "r" );
1937      for (int i=0; i<48; i++) {
1938        int   tid=0;
1939        float  tp[4];
1940        if(fscanf(fattin,"%d %f %f %f %f",
1941                  &tid, &tp[0], &tp[1], &tp[2], &tp[3])!=5) break;
1942        parPos[i].Set(4,tp);
1943      }
1944      fclose(fattin);
1945    }
1946    
1947    //------------------------------------------------------------------------
1948    void ToFdEdx::ReadParAtt( const char *fname )
1949    {
1950      printf("read %s\n",fname);
1951      FILE *fattin = fopen( fname , "r" );
1952      for (int i=0; i<48; i++) {
1953        int   tid=0;
1954        float  tp[6];
1955        if(fscanf(fattin,"%d %f %f %f %f %f %f",
1956                  &tid, &tp[0], &tp[1], &tp[2], &tp[3], &tp[4], &tp[5] )!=7) break;
1957        parAtt[i].Set(6,tp);
1958      }
1959      fclose(fattin);
1960    }
1961    
1962    
1963    
1964    
1965    
1966    
1967    double ToFdEdx::f_att( TArrayF &p, float x )
1968    {
1969      return
1970        p[0] +
1971        p[1]*x +
1972        p[2]*x*x +
1973        p[3]*x*x*x +
1974        p[4]*x*x*x*x +
1975        p[5]*x*x*x*x*x;
1976    }
1977    //------------------------------------------------------------------------
1978    double ToFdEdx::f_att5B( float x )
1979    {
1980      return
1981        101.9409 +
1982        6.643781*x +
1983        0.2765518*x*x +
1984        0.004617647*x*x*x +
1985        0.0006195132*x*x*x*x +
1986        0.00002813734*x*x*x*x*x;
1987    }
1988    
1989    
1990    double ToFdEdx::f_pos( TArrayF &p, float x )
1991    {
1992      return
1993        p[0] +
1994        p[1]*x +
1995        p[2]*x*x +
1996        p[3]*x*x*x;
1997    }
1998    
1999    double ToFdEdx::f_pos5B( float x )
2000    {
2001      return
2002        15.45132 +
2003        0.8369721*x +
2004        0.0005*x*x;
2005    }
2006    
2007    
2008    
2009    double ToFdEdx::f_adcPC( float x )
2010    {
2011      return 28.12+0.6312*x-5.647e-05*x*x+3.064e-08*x*x*x;
2012    }
2013    
2014    
2015    float ToFdEdx::Get_adc_he( int id, float pl_x[6], float pl_y[6])
2016    {
2017    
2018      //
2019      // input: id - pmt [0:47}
2020      //             pl_x - coord x of the tof plane
2021      //             pl_y - coord y
2022    
2023       adc_he = 0;
2024      if( eGeom.GetXY(id)==1 )  adc_he = f_att( (parAtt[id]), pl_x[eGeom.GetPlane(id)] );
2025      if( eGeom.GetXY(id)==2 )  adc_he = f_att( (parAtt[id]), pl_y[eGeom.GetPlane(id)] );
2026      return adc_he;
2027    }
2028    
2029    //------------------------------------------------------------------------
2030    double ToFdEdx::f_BB( TArrayF &p, float x )
2031    {
2032      return  p[0]/(x*x)*(log(x*x/(1-x*x)) - p[1]*x*x - p[2]);
2033    }
2034    
2035    //------------------------------------------------------------------------
2036    double ToFdEdx::f_BB5B( float x )
2037    {
2038      return  0.165797/(x*x)*(log(x*x/(1-x*x)) + 140.481*x*x + 52.9258);
2039    }
2040    //------------------------------------------------------------------------
2041    double ToFdEdx::f_desatBB( TArrayF &p, float x )
2042    {
2043      return
2044        p[0] +
2045        p[1]*x +
2046        p[2]*x*x;
2047    }
2048    
2049    //------------------------------------------------------------------------
2050    double ToFdEdx::f_desatBB5B( float x )
2051    {
2052      return
2053        -2.4 +
2054        0.75*x +
2055        0.009*x*x;
2056    }
2057    
2058    
2059    
2060    
2061    
2062    
2063    
2064    
2065    
2066    
2067    
2068    
2069    
2070    
2071    
2072    
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.23