/[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.16 by mocchiut, Mon Apr 30 15:46:30 2007 UTC revision 1.33 by mocchiut, Tue Nov 9 15:21:02 2010 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 15  ToFPMT::ToFPMT(){ Line 20  ToFPMT::ToFPMT(){
20    pmt_id = 0;    pmt_id = 0;
21    adc = 0.;    adc = 0.;
22    tdc_tw = 0.;    tdc_tw = 0.;
23      tdc = 0.;
24  }  }
25    
26  ToFPMT::ToFPMT(const ToFPMT &t){  ToFPMT::ToFPMT(const ToFPMT &t){
27    pmt_id = t.pmt_id;    pmt_id = t.pmt_id;
28    adc = t.adc;    adc = t.adc;
29    tdc_tw = t.tdc_tw;    tdc_tw = t.tdc_tw;
30      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.;
37      tdc = 0.;
38  }  }
39    
40    
# Line 50  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 104  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 113  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 213  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 378  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 391  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 424  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 436  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 505  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 538  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 556  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 574  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 591  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 608  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 625  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 693  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 927  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;  
857    
858      int somma=0;
859    if(plane == 0){    int np=plane;
860      padid=paddle;    for(Int_t j=0; j<np; j++){
861    }      somma+=pads[j];
   
   if(plane == 1){  
     padid=pads11+paddle;  
862    }    }
863      padid=paddle+somma;
   if(plane == 2){  
     padid=pads11+pads12+paddle;  
   }  
   
   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 991  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 1046  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 1096  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      printf("\n\n\n ERROR: NOT IMPLEMENTED ANYMORE, write Emiliano if you need this method (Emiliano.Mocchiutti@ts.infn.it) \n\n\n");
1153      return(-1);
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    //  return(0);
1417    }
1418    
1419    
1420    ToFdEdx::ToFdEdx()
1421    {
1422      memset(conn,0,12*sizeof(Bool_t));
1423      memset(ts,0,12*sizeof(UInt_t));
1424      memset(te,0,12*sizeof(UInt_t));
1425      Define_PMTsat();
1426      Clear();
1427    }
1428    //------------------------------------------------------------------------
1429    void ToFdEdx::CheckConnectors(UInt_t atime, GL_PARAM *glparam, TSQLServer *dbc)
1430    {
1431      for(int i=0; i<12; i++){
1432        if(atime<=ts[i] || atime>te[i]){
1433          Int_t error=glparam->Query_GL_PARAM(atime,210+i,dbc); // parameters stored in DB in GL_PRAM table
1434          if ( error<0 ) {
1435            conn[i]=false;
1436            ts[i]=0;
1437            te[i]=numeric_limits<UInt_t>::max();
1438          };
1439          if ( !error ){
1440            conn[i]=true;
1441            ts[i]=glparam->FROM_TIME;
1442            te[i]=glparam->TO_TIME;
1443          }
1444          if ( error>0 ){
1445            conn[i]=false;
1446            ts[i]=glparam->TO_TIME;
1447            TSQLResult *pResult;
1448            TSQLRow *row;
1449            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);
1450            pResult=dbc->Query(query.Data());
1451            if(!pResult->GetRowCount()){
1452              te[i]=numeric_limits<UInt_t>::max();
1453            }else{
1454              row=pResult->Next();
1455              te[i]=(UInt_t)atoll(row->GetField(0));
1456            }
1457          }
1458          //
1459          
1460        }
1461      }
1462    
1463    }
1464    //------------------------------------------------------------------------
1465    void ToFdEdx::Clear(Option_t *option)
1466    {
1467      //
1468      // Set arrays and initialize structure
1469      eDEDXpmt.Set(48);    eDEDXpmt.Reset(-1);   // Set array size  and reset structure
1470      //
1471    };
1472    
1473    //------------------------------------------------------------------------
1474    void ToFdEdx::Print(Option_t *option)
1475    {
1476      //
1477      printf("========================================================================\n");
1478    
1479    };
1480    
1481    //------------------------------------------------------------------------
1482    void ToFdEdx::Init(pamela::tof::TofEvent *tofl0)
1483    {
1484      //
1485      ToFLevel2 tf;
1486      for (Int_t gg=0; gg<4;gg++){
1487        for (Int_t hh=0; hh<12;hh++){
1488          //          tofinput_.tdc[hh][gg]=tofEvent->tdc[gg][hh];          
1489          int mm = tf.GetPMTid(gg,hh);        
1490          adc[mm]=tofl0->adc[gg][hh];
1491        };      
1492      };
1493      
1494    };
1495    
1496    //------------------------------------------------------------------------
1497    void ToFdEdx::Init(Int_t gg, Int_t hh, Float_t adce)
1498    {
1499      //
1500      ToFLevel2 tf;
1501      //  for (Int_t gg=0; gg<4;gg++){
1502      //    for (Int_t hh=0; hh<12;hh++){
1503      int mm = tf.GetPMTid(gg,hh);    
1504      adc[mm]=adce;
1505      
1506    };
1507    //------------------------------------------------------------------------
1508    void ToFdEdx::Process(UInt_t atime, Float_t betamean, Float_t *xtr_tof, Float_t *ytr_tof, Int_t exitat)
1509    {
1510      // the parameters should be already initialised by InitPar()
1511      Clear();
1512    
1513     // define angle:  
1514      double dx   = xtr_tof[1] - xtr_tof[5];
1515      double dy   = ytr_tof[0] - ytr_tof[4];
1516      double dr   = sqrt(dx*dx+dy*dy);
1517      double theta=atan(dr/76.81);
1518      //
1519      if ( xtr_tof[1] > 99. ||  xtr_tof[5] > 99. || ytr_tof[0] > 99. ||  ytr_tof[4] > 99. ) theta = 0.;
1520      for (Int_t ii=0; ii<6; ii++){
1521        if ( xtr_tof[ii] > 99. ) xtr_tof[ii] = 0.;
1522        if ( ytr_tof[ii] > 99. ) ytr_tof[ii] = 0.;
1523      };
1524      //
1525      
1526    
1527      //--------------------- TABLE OF PERIODS WITH HV PROBLEMS ----------------------------
1528      
1529      int Aconn=conn[0];    // PMT 0,20,22,24
1530      int Bconn=conn[1];    // PMT 6,12,26,34
1531      int Cconn=conn[2];    // PMT 4,14,28,32
1532      int Dconn=conn[3];    // PMT 2,8,10,30
1533      int Econn=conn[4];    // PMT 42,43,44,47
1534      int Fconn=conn[5];    // PMT 7,19,23,27
1535      int Gconn=conn[6];    // PMT 3,11,25,33
1536      int Hconn=conn[7];    // PMT 1,9,13,21
1537      int Iconn=conn[8];    // PMT 5,29,31,35
1538      int Lconn=conn[9];    // PMT 37,40,45,46
1539      int Mconn=conn[10];    // PMT 15,16,17,18
1540      int Nconn=conn[11];    // PMT 36,38,39,41
1541      if( false ) cout << Gconn << Iconn << Lconn <<endl; // to avoid compilation warnings
1542        
1543      for( int ii=0; ii<48; ii++ ) {
1544        //
1545        //    printf(" ii %i beta %f atime %u xtr 1 %f ytr 1 %f adc %f \n",ii,betamean,atime,xtr_tof[0],ytr_tof[0],adc[ii]);
1546    
1547        if( adc[ii] >= 4095 ){
1548          eDEDXpmt[ii] = 0.;
1549          continue; // EMILIANO
1550        };
1551    
1552        if( adc[ii] >= (PMTsat[ii]-5) && adc[ii] != 4095 ){
1553          eDEDXpmt[ii] = 1000.;
1554          continue; // EMILIANO
1555        };
1556    
1557        if( adc[ii] <= 0. ) {
1558          eDEDXpmt[ii] = 1500.;
1559          continue;
1560        };
1561        //
1562        double adcpC   = f_adcPC( adc[ii] );    // - adc conversion in pC
1563        if ( exitat == 0 ){
1564          eDEDXpmt[ii]=(Float_t)adcpC;
1565          continue;
1566        }
1567    
1568        double adccorr = adcpC*fabs(cos(theta));    
1569        if(adccorr<=0.)           continue;
1570        if ( exitat == 1 ){
1571          eDEDXpmt[ii]=(Float_t)adccorr;
1572          continue;
1573        }
1574    
1575        //    int standard=0;
1576        int S115B_ok=0;
1577        int S115B_break=0;
1578    
1579        if(atime<1158720000)S115B_ok=1;
1580        else S115B_break=1;
1581    
1582    
1583        //------------------------------------------------------------------------
1584    
1585        //---------------------------------------------------- Z reconstruction
1586    
1587        double adcHe, adcnorm, adclin, dEdx, Zeta;
1588    
1589        adcHe=-2;
1590        adcnorm=-2;
1591        adclin=-2;
1592        dEdx=-2;
1593        Zeta=-2;
1594        Double_t correction = 1.;
1595    
1596        if(Aconn==1 && (ii==0 || ii==20 || ii==22 || ii==24)){
1597          correction = 1.675;
1598        }
1599        else if(Bconn==1 && (ii==6 || ii==12 || ii==26 || ii==34)){
1600          correction = 2.482;
1601        }
1602        else if(Cconn==1 && (ii==4 || ii==14 || ii==28 || ii==32)){
1603          correction = 1.464;
1604        }
1605        else if(Dconn==1 && (ii==2 || ii==8 || ii==10 || ii==30)){
1606          correction = 1.995;
1607        }
1608        else if(Econn==1 && (ii==42 || ii==43 || ii==44 || ii==47)){
1609          correction = 1.273;
1610        }
1611        else if(Fconn==1 && (ii==7 || ii==19 || ii==23 || ii==27)){
1612          correction = 1.565;
1613        }
1614        else if(Mconn==1 && (ii==15 || ii==16 || ii==17 || ii==18)){
1615          correction = 1.565;
1616        }
1617        else if(Nconn==1 && (ii==36 || ii==38 || ii==39 || ii==41)){
1618          correction = 1.018;
1619        }
1620        else if(Hconn==1 && (ii==1 || ii==13 || ii==21 || (ii==9&&S115B_ok==1))){
1621          correction = 1.84;
1622        }
1623        else if(S115B_break==1 && ii==9 && Hconn==1){
1624          correction = 1.64;
1625        }
1626        else correction = 1.;
1627        
1628        if( S115B_break==1 ){
1629          adcHe   = f_att5B( ytr_tof[0] )/correction;
1630        } else {
1631          adcHe   = Get_adc_he(ii, xtr_tof, ytr_tof)/correction;
1632        };
1633        if(adcHe<=0)   continue;
1634        if ( exitat == 2 ){
1635          if(ii==9 && S115B_break==1)  eDEDXpmt[ii]=36.*(Float_t)adccorr/adcHe;
1636          else  adclin  = 4.*(Float_t)adccorr/adcHe;
1637          continue;
1638        }
1639    
1640        if(ii==9 && S115B_break==1)  adcnorm = f_pos5B(adccorr);
1641        else adcnorm = f_pos( (parPos[ii]), adccorr);
1642        if(adcnorm<=0) continue;
1643        if(ii==9 && S115B_break==1)  adclin  = 36.*adcnorm/adcHe;
1644        else  adclin  = 4.*adcnorm/adcHe;
1645        if(adclin<=0)  continue;
1646        if ( exitat == 3 ){
1647          if(ii==9 && S115B_break==1)  eDEDXpmt[ii]=(Float_t)adclin;
1648          else  eDEDXpmt[ii]=(Float_t)adclin;
1649          continue;
1650        }
1651        //
1652        if ( betamean > 99. ){
1653          //      eDEDXpmt.AddAt((Float_t)adclin,ii);
1654          eDEDXpmt[ii]=(Float_t)adclin;
1655          //      printf(" AAPMT IS %i dedx is %f vector is %f \n",ii,adclin,eDEDXpmt[ii]);
1656          continue;
1657        };
1658        //
1659        double dEdxHe=-2;
1660        if(ii==9 && S115B_break==1){
1661          if( betamean <1. ) dEdxHe = f_BB5B( betamean );
1662          else                       dEdxHe = 33;
1663        } else {
1664          if( betamean <1. ) dEdxHe = f_BB( (parBBneg[ii]), betamean );
1665          else                       dEdxHe = parBBpos[ii];
1666        }
1667        
1668        
1669        if(dEdxHe<=0){
1670          eDEDXpmt[ii]=(Float_t)adclin;
1671          continue;
1672        };
1673    
1674        if(ii==9 && S115B_break==1)  dEdx = f_desatBB5B( adclin );
1675        else  dEdx = f_desatBB((parDesatBB[ii]), adclin );
1676    
1677        if(dEdx<=0){
1678          eDEDXpmt[ii]=(Float_t)adclin;
1679          continue;
1680        };
1681    
1682        eDEDXpmt[ii]=(Float_t)dEdx;
1683        //    eDEDXpmt.AddAt((Float_t)dEdx,ii);
1684    
1685        //    printf(" PMT IS %i dedx is %f vector is %f \n",ii,dEdx,eDEDXpmt[ii]);
1686    
1687      }  //end loop on 48 PMT
1688    
1689    };
1690    
1691    
1692    //------------------------------------------------------------------------
1693    void ToFdEdx::Define_PMTsat()
1694    {
1695      Float_t  sat[48] = {
1696        3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
1697        3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
1698        3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
1699        3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
1700        3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
1701        3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };
1702      PMTsat.Set(48,sat);
1703    }
1704    
1705    //------------------------------------------------------------------------
1706    void ToFdEdx::ReadParBBpos( const char *fname )
1707    {
1708      //  printf("read %s\n",fname);
1709      parBBpos.Set(48);
1710      FILE *fattin = fopen( fname , "r" );
1711      for (int i=0; i<48; i++) {
1712        int   tid=0;
1713        float  tp;
1714        if(fscanf(fattin,"%d %f",
1715                  &tid, &tp )!=2) break;
1716        parBBpos[i]=tp;
1717      }
1718      fclose(fattin);
1719    }
1720    
1721    //------------------------------------------------------------------------
1722    void ToFdEdx::ReadParDesatBB( const char *fname )
1723    {
1724      //  printf("read %s\n",fname);
1725      FILE *fattin = fopen( fname , "r" );
1726      for (int i=0; i<48; i++) {
1727        int   tid=0;
1728        float  tp[3];
1729        if(fscanf(fattin,"%d %f %f %f",
1730                  &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
1731        parDesatBB[i].Set(3,tp);
1732      }
1733      fclose(fattin);
1734    }
1735    
1736    
1737    //------------------------------------------------------------------------
1738    void ToFdEdx::ReadParBBneg( const char *fname )
1739    
1740    {
1741      //  printf("read %s\n",fname);
1742      FILE *fattin = fopen( fname , "r" );
1743      for (int i=0; i<48; i++) {
1744        int   tid=0;
1745        float  tp[3];
1746        if(fscanf(fattin,"%d %f %f %f",
1747                  &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
1748        parBBneg[i].Set(3,tp);
1749      }
1750      fclose(fattin);
1751    }
1752    
1753    //------------------------------------------------------------------------
1754    void ToFdEdx::ReadParPos( const char *fname )
1755    {
1756      //  printf("read %s\n",fname);
1757      FILE *fattin = fopen( fname , "r" );
1758      for (int i=0; i<48; i++) {
1759        int   tid=0;
1760        float  tp[4];
1761        if(fscanf(fattin,"%d %f %f %f %f",
1762                  &tid, &tp[0], &tp[1], &tp[2], &tp[3])!=5) break;
1763        parPos[i].Set(4,tp);
1764      }
1765      fclose(fattin);
1766    }
1767    
1768    //------------------------------------------------------------------------
1769    void ToFdEdx::ReadParAtt( const char *fname )
1770    {
1771      //  printf("read %s\n",fname);
1772      FILE *fattin = fopen( fname , "r" );
1773      for (int i=0; i<48; i++) {
1774        int   tid=0;
1775        float  tp[6];
1776        if(fscanf(fattin,"%d %f %f %f %f %f %f",
1777                  &tid, &tp[0], &tp[1], &tp[2], &tp[3], &tp[4], &tp[5] )!=7) break;
1778        parAtt[i].Set(6,tp);
1779      }
1780      fclose(fattin);
1781    }
1782    
1783    
1784    
1785    
1786    
1787    
1788    double ToFdEdx::f_att( TArrayF &p, float x )
1789    {
1790      return
1791        p[0] +
1792        p[1]*x +
1793        p[2]*x*x +
1794        p[3]*x*x*x +
1795        p[4]*x*x*x*x +
1796        p[5]*x*x*x*x*x;
1797    }
1798    //------------------------------------------------------------------------
1799    double ToFdEdx::f_att5B( float x )
1800    {
1801      return
1802        101.9409 +
1803        6.643781*x +
1804        0.2765518*x*x +
1805        0.004617647*x*x*x +
1806        0.0006195132*x*x*x*x +
1807        0.00002813734*x*x*x*x*x;
1808    }
1809    
1810    
1811    double ToFdEdx::f_pos( TArrayF &p, float x )
1812    {
1813      return
1814        p[0] +
1815        p[1]*x +
1816        p[2]*x*x +
1817        p[3]*x*x*x;
1818    }
1819    
1820    double ToFdEdx::f_pos5B( float x )
1821    {
1822      return
1823        15.45132 +
1824        0.8369721*x +
1825        0.0005*x*x;
1826    }
1827    
1828    
1829    
1830    double ToFdEdx::f_adcPC( float x )
1831    {
1832      return 28.12+0.6312*x-5.647e-05*x*x+3.064e-08*x*x*x;
1833    }
1834    
1835    
1836    float ToFdEdx::Get_adc_he( int id, float pl_x[6], float pl_y[6])
1837    {
1838    
1839      //
1840      // input: id - pmt [0:47}
1841      //             pl_x - coord x of the tof plane
1842      //             pl_y - coord y
1843    
1844      adc_he = 0;
1845      if( eGeom.GetXY(id)==1 )  adc_he = f_att( (parAtt[id]), pl_x[eGeom.GetPlane(id)] );
1846      if( eGeom.GetXY(id)==2 )  adc_he = f_att( (parAtt[id]), pl_y[eGeom.GetPlane(id)] );
1847      return adc_he;
1848    }
1849    
1850    //------------------------------------------------------------------------
1851    double ToFdEdx::f_BB( TArrayF &p, float x )
1852    {
1853      return  p[0]/(x*x)*(log(x*x/(1-x*x)) - p[1]*x*x - p[2]);
1854    }
1855    
1856    //------------------------------------------------------------------------
1857    double ToFdEdx::f_BB5B( float x )
1858    {
1859      return  0.165797/(x*x)*(log(x*x/(1-x*x)) + 140.481*x*x + 52.9258);
1860    }
1861    //------------------------------------------------------------------------
1862    double ToFdEdx::f_desatBB( TArrayF &p, float x )
1863    {
1864      return
1865        p[0] +
1866        p[1]*x +
1867        p[2]*x*x;
1868    }
1869    
1870    //------------------------------------------------------------------------
1871    double ToFdEdx::f_desatBB5B( float x )
1872    {
1873      return
1874        -2.4 +
1875        0.75*x +
1876        0.009*x*x;
1877    }
1878    
1879    
1880    
1881    
1882    
1883    
1884    
1885    
1886    
1887    
1888    
1889    
1890    
1891    
1892    
1893    
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.23