/[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.1 by mocchiut, Sat Jun 17 12:14:56 2006 UTC revision 1.47 by pam-fi, Tue May 19 13:46:46 2015 UTC
# Line 1  Line 1 
1  #include <TObject.h>  /**
2     * \file ToFLevel2.cpp
3     * \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    
11  #include <ToFLevel2.h>  #include <ToFLevel2.h>
 #include <iostream>  
12  using namespace std;  using namespace std;
13    ClassImp(ToFPMT);
14    ClassImp(ToFdEdx);
15    ClassImp(ToFGeom);
16  ClassImp(ToFTrkVar);  ClassImp(ToFTrkVar);
17  ClassImp(ToFLevel2);  ClassImp(ToFLevel2);
18    
19  ToFTrkVar::ToFTrkVar() {  ToFPMT::ToFPMT(){
20      pmt_id = 0;
21    trkseqno = 0;    adc = 0.;
22      tdc_tw = 0.;
23      tdc = 0.;
24      l0flag_adc = 0.;
25      l0flag_tdc = 0.;
26    }
27    
28    for (Int_t kk=0; kk<13;kk++){  ToFPMT::ToFPMT(const ToFPMT &t){
29      beta_a[kk] = 0;    pmt_id = t.pmt_id;
30    }    adc = t.adc;
31      tdc_tw = t.tdc_tw;
32      tdc = t.tdc;
33    }
34    
35    for (Int_t kk=0; kk<4;kk++){  void ToFPMT::Clear(Option_t *t){
36      for (Int_t hh=0; hh<12;hh++){    pmt_id = 0;
37        adc_c[hh][kk] = 0;    adc = 0.;
38      }    tdc_tw = 0.;
39    }    tdc = 0.;
40  }  }
41    
 ToFTrkVar::ToFTrkVar(const ToFTrkVar &t){  
42    
   trkseqno = t.trkseqno;    
43    
44    memcpy(adc_c,t.adc_c,sizeof(adc_c));  ToFTrkVar::ToFTrkVar() {
45    memcpy(beta_a,t.beta_a,sizeof(beta_a));    trkseqno = 0;
46  }    npmttdc = 0;
47      npmtadc = 0;
48      pmttdc = TArrayI(48);
49      pmtadc = TArrayI(48);
50      tdcflag = TArrayI(48); // gf: 30 Nov 2006
51      adcflag = TArrayI(48); // gf: 30 Nov 2006
52      dedx = TArrayF(48);
53      //
54      //
55      memset(beta,  0, 13*sizeof(Float_t));
56      memset(xtofpos,  0, 3*sizeof(Float_t));
57      memset(ytofpos,  0, 3*sizeof(Float_t));
58      memset(xtr_tof,  0, 6*sizeof(Float_t));
59      memset(ytr_tof,  0, 6*sizeof(Float_t));
60      //
61    };
62    
63  ToFLevel2::ToFLevel2() {      void ToFTrkVar::Clear(Option_t *t) {
64      trkseqno = 0;
65      npmttdc = 0;
66      npmtadc = 0;
67      pmttdc.Reset();
68      pmtadc.Reset();
69      tdcflag.Reset(); // gf: 30 Nov 2006
70      adcflag.Reset(); // gf: 30 Nov 2006
71      dedx.Reset();
72    //    //
73    ToFTrk = new TClonesArray("ToFTrkVar",1);    memset(beta,  0, 13*sizeof(Float_t));
74      memset(xtofpos,  0, 3*sizeof(Float_t));
75      memset(ytofpos,  0, 3*sizeof(Float_t));
76      memset(xtr_tof,  0, 6*sizeof(Float_t));
77      memset(ytr_tof,  0, 6*sizeof(Float_t));
78    //    //
79    for (Int_t kk=0; kk<3;kk++){  };
     xtofpos[kk] = 0.;  
     ytofpos[kk] = 0.;  
   }  
80    
81    for (Int_t kk=0; kk<6;kk++){  ToFTrkVar::ToFTrkVar(const ToFTrkVar &t){
     tof_i_flag[kk] = 0;  
     tof_j_flag[kk] = 0;  
   }  
     
   for (Int_t kk=0; kk<13;kk++){  
     betatof_a[kk] = 0;  
   }  
82    
83    for (Int_t kk=0; kk<4;kk++){    trkseqno = t.trkseqno;  
84      for (Int_t hh=0; hh<12;hh++){    //
85        adctof_c[hh][kk] = 0;    npmttdc = t.npmttdc;
86      }    npmtadc = t.npmtadc;
87    }      (t.pmttdc).Copy(pmttdc);
88      (t.pmtadc).Copy(pmtadc);
89      (t.tdcflag).Copy(tdcflag); // gf: 30 Nov 2006
90      (t.adcflag).Copy(adcflag); // gf: 30 Nov 2006
91      (t.dedx).Copy(dedx);
92      //
93      memcpy(beta,t.beta,sizeof(beta));
94      memcpy(xtofpos,t.xtofpos,sizeof(xtofpos));
95      memcpy(ytofpos,t.ytofpos,sizeof(ytofpos));
96      memcpy(xtr_tof,t.xtr_tof,sizeof(xtr_tof));
97      memcpy(ytr_tof,t.ytr_tof,sizeof(ytr_tof));
98      //
99    };
100    
101    for (Int_t kk=0; kk<4;kk++){  ToFLevel2::ToFLevel2() {    
102      for (Int_t hh=0; hh<12;hh++){    //
103        tdc_c[hh][kk] = 0;  //  PMT = new TClonesArray("ToFPMT",12); //ELENA
104      }  //  ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
105    }      PMT = 0; //ELENA
106      ToFTrk = 0; //ELENA
107      //
108      this->Clear();
109      //
110    };
111    
112    void ToFLevel2::Set(){//ELENA
113        if(!PMT)PMT = new TClonesArray("ToFPMT",12); //ELENA
114        if(!ToFTrk)ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
115    }//ELENA
116    //--------------------------------------
117    //
118    //
119    //--------------------------------------
120    void ToFLevel2::SetTrackArray(TClonesArray *track){//ELENA
121        if(track && strcmp(track->GetClass()->GetName(),"ToFTrkVar")==0){
122            if(ToFTrk)ToFTrk->Clear("C");
123            ToFTrk = track;
124        }
125    }
126    
127    void ToFLevel2::Clear(Option_t *t){
128      //
129      if(ToFTrk)ToFTrk->Delete(); //ELENA
130      if(PMT)PMT->Delete(); //ELENA
131      memset(tof_j_flag, 0, 6*sizeof(Int_t));
132      unpackError = 0;
133      unpackWarning = 0;
134      //
135  };  };
136    
137    void ToFLevel2::Delete(Option_t *t){ //ELENA
138      //
139      if(ToFTrk){
140          ToFTrk->Delete(); //ELENA
141          delete ToFTrk;  //ELENA
142      }
143      if(PMT){
144          PMT->Delete(); //ELENA
145          delete PMT; //ELENA
146      } //ELENA
147      //
148    }; //ELENA
149    
150    /**
151     * Retrieves the itrk-th tof track stored in the array
152     * @param itrk Array index (min 0, max ToFLevel2::ntrk())
153     *
154     */
155  ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t itrk){  ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t itrk){
156    //        //    
157    if(itrk >= ntrk()){    if(itrk >= ntrk()){
# Line 69  ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t Line 160  ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t
160      return(NULL);      return(NULL);
161    }      }  
162    //    //
163      if(!ToFTrk)return 0; //ELENA
164    TClonesArray &t = *(ToFTrk);    TClonesArray &t = *(ToFTrk);
165    ToFTrkVar *toftrack = (ToFTrkVar*)t[itrk];    ToFTrkVar *toftrack = (ToFTrkVar*)t[itrk];
166    return toftrack;    return toftrack;
167  }  }
168    
169    /**
170     * Retrieves the tof track matching the seqno-th tracker stored track.
171     * @param seqno Track sequential number
172     * (seqno = -1 for standalone tof track, seqno=0-TrkLevel2::ntrk() for tof tracks associated to a tracker track)
173     *
174     */
175    ToFTrkVar *ToFLevel2::GetToFStoredTrack(int seqno){
176    
177      if( ntrk()==0 ){
178        printf("ToFLevel2::GetToFStoredTrack(int) : requested tracker SeqNo %i but no ToFrimeter tracks are stored\n",seqno);
179        return NULL;
180      };
181      
182      ToFTrkVar *c = 0;
183      Int_t it_tof=0;
184        
185      do {
186        c = GetToFTrkVar(it_tof);
187        it_tof++;
188      } while( c && seqno != c->trkseqno && it_tof < ntrk());      
189      
190      if(!c || seqno != c->trkseqno){
191        c = 0;
192        if(seqno!=-1 ) printf("ToFLevel2::GetToFStoredTrack(int) : requested tracker SeqNo %i does not match ToFrimeter stored tracks\n",seqno);
193      };
194      return c;
195        
196    }
197    
198    
199    ToFPMT *ToFLevel2::GetToFPMT(Int_t ihit){
200      //    
201      if(ihit >= npmt()){
202        printf(" ToFLevel2 ERROR: pmt variables set %i does not exists! \n",ihit);
203        printf("                  stored pmt variables = %i \n",npmt());
204        return(NULL);
205      }  
206      //
207      if(!PMT)return 0; //ELENA
208      TClonesArray &t = *(PMT);
209      ToFPMT *tofpmt = (ToFPMT*)t[ihit];
210      return tofpmt;
211    }
212  //--------------------------------------  //--------------------------------------
213  //  //
214  //  //
215  //--------------------------------------  //--------------------------------------
216  /**  /**
217   * Method to get the paddle ID (11 12 21 22 31 32) from the paddle index (0 1 2 3 4 5)   * Method to get the plane ID (11 12 21 22 31 32) from the plane index (0 1 2 3 4 5)
218     * @param Plane index (0,1,2,3,4,5).
219   */   */
220    Int_t  ToFLevel2::GetToFPlaneID(Int_t ip){    Int_t  ToFLevel2::GetToFPlaneID(Int_t ip){
221        if(ip>=0 && ip<6)return 10*((int)(ip/2+1.1))+(ip%2)+1;        if(ip>=0 && ip<6)return 10*((int)(ip/2+1.1))+(ip%2)+1;
222        else return -1;        else return -1;
223    };    };
224  /**  /**
225   * Method to get the paddle index (0 1 2 3 4 5) from the paddle ID (11 12 21 22 31 32)   * Method to get the plane index (0 1 2 3 4 5) from the plane ID (11 12 21 22 31 32)
226     * @param plane Plane ID (11, 12, 21, 22, 31, 32)
227   */   */
228    Int_t  ToFLevel2::GetToFPlaneIndex(Int_t plane_id){    Int_t  ToFLevel2::GetToFPlaneIndex(Int_t plane_id){
229        if(        if(
# Line 99  ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t Line 237  ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t
237        else return -1;        else return -1;
238    };    };
239  /**  /**
240   * Method to know if a given ToF paddle was hit, that is there is a TDC signal from both PMTs   * Method to know if a given ToF paddle was hit, that is there is a TDC signal
241     * from both PMTs. The method uses the "tof_j_flag" variable.
242   * @param plane Plane ID (11, 12, 21, 22, 31, 32) or Plane index (0,1,2,3,4,5).   * @param plane Plane ID (11, 12, 21, 22, 31, 32) or Plane index (0,1,2,3,4,5).
243   * @param paddle_id Paddle ID.   * @param paddle_id Paddle ID.
244   * @return 1 if the paddle was hit.   * @return 1 if the paddle was hit.
# Line 109  ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t Line 248  ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t
248      if     (plane>=6             ) ip = GetToFPlaneIndex(plane);      if     (plane>=6             ) ip = GetToFPlaneIndex(plane);
249      else if(plane>=0 && plane < 6) ip = plane;      else if(plane>=0 && plane < 6) ip = plane;
250      Int_t flag=0;      Int_t flag=0;
251      if(ip != -1)flag = tof_j_flag[ip] & (int)pow(2,paddle_id);      if(ip != -1)flag = tof_j_flag[ip] & (int)pow(2.,(double)paddle_id);
252      if(      if(
253         (ip == 0 && paddle_id < 8 && flag) ||         (ip == 0 && paddle_id < 8 && flag) ||
254         (ip == 1 && paddle_id < 6 && flag) ||         (ip == 1 && paddle_id < 6 && flag) ||
# Line 120  ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t Line 259  ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t
259         false) return true;         false) return true;
260      else return false;      else return false;
261  };  };
262    
263  /**  /**
264   * Method to get the number of hit paddles on a ToF plane.   * Strict method to get the number of hit paddles on a ToF plane.
265     * The method uses "HitPaddle" which checks if there is a TDC signal
266     * from both PMTs.
267   * @param plane Plane ID (11, 12, 21, 22, 31, 32) or Plane index (0,1,2,3,4,5).   * @param plane Plane ID (11, 12, 21, 22, 31, 32) or Plane index (0,1,2,3,4,5).
268   */   */
269  Int_t ToFLevel2::GetNHitPaddles(Int_t plane){  Int_t ToFLevel2::GetNHitPaddles(Int_t plane){
# Line 129  Int_t ToFLevel2::GetNHitPaddles(Int_t pl Line 271  Int_t ToFLevel2::GetNHitPaddles(Int_t pl
271      for(Int_t i=0; i<8; i++)npad = npad + (int)HitPaddle(plane,i);      for(Int_t i=0; i<8; i++)npad = npad + (int)HitPaddle(plane,i);
272      return npad;      return npad;
273  };  };
274    
275    /**
276     * Optional method to get the number of hit paddles on a ToF plane.
277     * The method does NOT check if there is a signal from both PMTs, it only
278     * checks if there is some PMT signal in a paddle
279     * @param plane Plane ID (11, 12, 21, 22, 31, 32) or Plane index (0,1,2,3,4,5).
280     */
281    Int_t ToFLevel2::GetTrueNHitPaddles(Int_t plane){
282        Int_t npad=0;
283        TClonesArray* Pmt = this->PMT;
284        int paddle[24];
285        memset(paddle,0, 24*sizeof(int));
286        for(int i=0; i<Pmt->GetEntries(); i++) {  //loop per vedere quale TOF è colpito
287          ToFPMT* pmthit = (ToFPMT*)Pmt->At(i);
288          int pplane = -1;
289          int ppaddle = -1;
290          GetPMTPaddle(pmthit->pmt_id,pplane,ppaddle);
291          if ( pplane == plane ) paddle[ppaddle]++;
292        }
293        for(int i=0;i<24;i++) if ( paddle[i]>0 ) npad++;
294    
295        return npad;
296    };
297    
298    //new, wm Feb 15
299    //wm Nov 08
300    //gf Apr 07
301    /**
302     * Method to get the mean dEdx from a ToF layer
303     * By definition there should be PMTs with dEdx values only in one paddle of a layer
304     * (the paddle hitted by the track), this method looks for the hitted paddle
305     * and gives the mean dEdx of that paddle as the output
306     * The method was modified for the "ToF-standalone" part in february 2015
307     * The "adcfl" option is not very useful (an artificial dEdx is per
308     * definition= 1 mip and not a real measurement), anyway left in the code
309     * @param notrack Track Number
310     * @param plane Plane index (0,1,2,3,4,5)
311     * @param adcflag in the plane (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )
312     */
313    Float_t ToFLevel2::GetdEdx(Int_t notrack, Int_t plane, Int_t adcfl){
314    //  ToFTrkVar *trk = GetToFTrkVar(notrack);
315        ToFTrkVar *trk = GetToFStoredTrack(notrack);//Elena 2015
316      return this->GetdEdx(trk, plane, adcfl);
317    }
318    
319    //new, wm Feb 15
320    //wm Nov 08
321    //gf Apr 07
322    /**
323     * Method to get the mean dEdx from a ToF layer
324     * By definition there should be PMTs with dEdx values only in one paddle of a layer
325     * (the paddle hitted by the track), this method looks for the hitted paddle
326     * and gives the mean dEdx of that paddle as the output
327     * The method was modified for the "ToF-standalone" part in february 2015
328     * The "adcfl" option is not very useful (an artificial dEdx is per
329     * definition= 1 mip and not a real measurement), anyway left in the code
330     * @param trk Pointer to TofTrkVar object
331     * @param plane Plane index (0,1,2,3,4,5)
332     * @param adcflag in the plane (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )
333     */
334    Float_t ToFLevel2::GetdEdx(ToFTrkVar *trk, Int_t plane, Int_t adcfl){
335      
336      Float_t dedx = 0.;
337      Float_t PadEdx =0.;
338      Int_t SatWarning;
339      Int_t pad=-1;
340      //
341      if(!trk) cout << "ToFLevel2::GetdEdx(...) ---> NULL ToFTrkVar obj "<<endl;
342      if(!trk) return 0; //ELENA
343      //
344      // ToF standalone part
345      //
346      if ( trk->trkseqno == -1 ){
347        
348        //    ToFTrkVar *t_tof = trk;
349        
350        // Find the hitted paddle  (two good TDC values) using the tof_j_flag (from tofl2com.for)
351        
352        Int_t Ipaddle=-1;
353        // if tof_j_flag == 0: no paddle was hitted. Otherwise decode tof_j_flag to get the paddle
354        if (this->tof_j_flag[plane] > 0)  Ipaddle = (Int_t)log2(this->tof_j_flag[plane]) ;
355        
356        Ipaddle =  (Int_t)log2(this->tof_j_flag[plane]) ;
357        
358        // Get the dEdx of this paddle using "GetdEdxPaddle"
359        if (Ipaddle>-1) {
360          Int_t pad = GetPaddleid(plane,Ipaddle);
361          GetdEdxPaddle(trk, pad, adcfl, PadEdx, SatWarning);
362          dedx = PadEdx;
363        }
364        
365        // If there was no correct hitted paddle, but there was one (and only one) paddle with some
366        // PMT entries in the PMT-class (found with "GetTrueNHitPaddles", use the dEdx of this paddle
367        
368        if ((Ipaddle<0) && (GetTrueNHitPaddles(plane)==1)) {
369          // find the paddle by looping over the paddles in each layer
370          // since GetTrueNHitPaddles==1 this is OK
371          for (Int_t ii=0; ii<GetNPaddle(plane); ii++){
372            Int_t paddleid=ii;
373            Int_t pad = GetPaddleid(plane,paddleid);
374            GetdEdxPaddle(trk, pad, adcfl, PadEdx, SatWarning);
375            dedx += PadEdx;
376          }
377        }
378      } else {
379        // track dependent dEdx: simple, there will be only one paddle hitted in    each layer
380        // so just loop over the paddles in each layer
381        for (Int_t ii=0; ii<GetNPaddle(plane); ii++){
382          Int_t paddleid=ii;
383          pad = GetPaddleid(plane,paddleid);
384          GetdEdxPaddle(trk, pad, adcfl, PadEdx, SatWarning);
385          dedx += PadEdx;
386        }
387      }
388      //
389      return(dedx);
390    }
391    
392    /**
393     * Method to fill the ADC_C 4x12 matrix with the dEdx values and the TDC 4x12 matrix
394     * with the time-walk corrected TDC values.
395     * @param notrack Track Number
396     * @param adc  ADC_C matrix with dEdx values
397     * @param tdc  TDC matrix
398     */
399    void ToFLevel2::GetMatrix(Int_t notrack, Float_t adc[4][12], Float_t tdc[4][12]){
400      //
401      for (Int_t aa=0; aa<4;aa++){
402        for (Int_t bb=0; bb<12;bb++){
403          adc[aa][bb] = 1000.;
404          tdc[aa][bb] = 4095.;
405        };
406      };
407      //
408      Int_t pmt_id = 0;
409      Int_t hh = 0;
410      Int_t kk = 0;
411      //
412      ToFTrkVar *trk = GetToFTrkVar(notrack);
413      if(!trk)return; //ELENA
414      //
415      for (Int_t i=0; i<trk->npmtadc; i++){
416        //
417        pmt_id = (trk->pmtadc).At(i);
418        //
419        GetPMTIndex(pmt_id,hh,kk);
420        adc[kk][hh] = (trk->dedx).At(i);  
421        //
422      };
423      //
424      for (Int_t i=0; i<npmt(); i++){
425        //
426        ToFPMT *pmt = GetToFPMT(i);
427        if(!pmt)break; //ELENA
428        //
429        GetPMTIndex(pmt->pmt_id,hh,kk);
430        //
431        tdc[kk][hh] = pmt->tdc_tw;  
432        //
433      };
434      //
435      return;
436    };
437    
438    
439    /**
440     * Method to get the plane index (0 - 5) for the PMT_ID as input
441     * @param pmt_id  PMT_ID (0 - 47)
442     */
443    Int_t ToFLevel2::GetPlaneIndex(Int_t pmt_id){
444      TString pmtname = GetPMTName(pmt_id);
445      pmtname.Resize(3);
446      if ( !strcmp(pmtname,"S11") ) return(0);
447      if ( !strcmp(pmtname,"S12") ) return(1);
448      if ( !strcmp(pmtname,"S21") ) return(2);
449      if ( !strcmp(pmtname,"S22") ) return(3);
450      if ( !strcmp(pmtname,"S31") ) return(4);
451      if ( !strcmp(pmtname,"S32") ) return(5);
452      return(-1);
453    };
454    
455    
456    /**
457     * Method to get the PMT_ID if the index (4,12) is given. We have 4 channels on
458     * each of the 12 half-boards, this method decodes which PMT is cables to which
459     * channel.
460     * @param hh Channel
461     * @param kk HalfBoard
462     */
463    Int_t ToFLevel2::GetPMTid(Int_t hh, Int_t kk){
464      //
465      static const short tof[4][24] = {
466        {4, 4,  4,  4,  1,  1, 2, 2,  3,  3, 3, 3,  3,  3, 1, 1,  1,  1, 2, 3,  3, 3, 3,  4},
467        {1, 3,  5,  7, 10, 12, 2, 4,  2,  4, 6, 8, 10, 12, 1, 5,  3,  9, 7, 9, 11, 1, 5,  9},
468        {2, 2,  2,  2,  1,  1, 1, 1,  4,  4, 4, 4,  4,  4, 2, 1,  2,  1, 2, 2,  2, 3, 3,  4},
469        {6, 8, 12, 10,  8,  6, 4, 2, 12, 10, 8, 6,  4,  2, 9, 7, 11, 11, 5, 3,  1, 3, 7, 11}
470      };
471      //
472      Int_t ind = 0;
473      Int_t k = 0;
474      while (k < 24){
475        Int_t j = 0;
476        while (j < 2){
477          Int_t ch = tof[2*j][k]     - 1;
478          Int_t hb = tof[2*j + 1][k] - 1;      
479          /* tofEvent->tdc[ch][hb] */      
480          if( ch == hh && hb == kk ){      
481            ind = 2*k + j;
482            break;
483          };
484          j++;
485        };
486        k++;
487      };
488      return ind;
489    };
490    
491    
492    /**
493     * Method to get the PMT index if the PMT ID is given. This method is the
494     * "reverse" of method "GetPMTid"
495     * @param ind  PMT_ID (0 - 47)
496     * @param hb   HalfBoard
497     * @param ch   Channel
498     */
499    void ToFLevel2::GetPMTIndex(Int_t ind, Int_t &hb, Int_t &ch){
500      //
501      static const short tof[4][24] = {
502        {4, 4,  4,  4,  1,  1, 2, 2,  3,  3, 3, 3,  3,  3, 1, 1,  1,  1, 2, 3,  3, 3, 3,  4},
503        {1, 3,  5,  7, 10, 12, 2, 4,  2,  4, 6, 8, 10, 12, 1, 5,  3,  9, 7, 9, 11, 1, 5,  9},
504        {2, 2,  2,  2,  1,  1, 1, 1,  4,  4, 4, 4,  4,  4, 2, 1,  2,  1, 2, 2,  2, 3, 3,  4},
505        {6, 8, 12, 10,  8,  6, 4, 2, 12, 10, 8, 6,  4,  2, 9, 7, 11, 11, 5, 3,  1, 3, 7, 11}
506      };
507      //
508      Int_t k = 0;
509      while (k < 24){
510        Int_t j = 0;
511        while (j < 2){
512          /* tofEvent->tdc[ch][hb] */            
513          if( ind == 2*k + j ){
514            ch = tof[2*j][k]     - 1;
515            hb = tof[2*j + 1][k] - 1;      
516            return;
517          };
518          j++;
519        };
520        k++;
521      };
522      return;
523    };
524    
525    
526    
527    //  wm Nov 08 revision - saturation values included
528    /// gf Apr 07
529    /**
530     * Method to get the dEdx from a given ToF paddle.
531     * If two PMTs are good, the mean dEdx of both PMTs is taken, otherwise
532     * just the dEdx of the "good" PMT. If both PMTs are above saturation => dEdx=1000
533     * @param notrack Track Number
534     * @param Paddle index (0,1,...,23).
535     * @param adcflag in the paddle (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )
536     * @param PadEdx dEdx from a given ToF paddle
537     * @param SatWarning 1 if the PMT ios near saturation region (adcraw ~3000)
538     */
539    void ToFLevel2::GetdEdxPaddle(Int_t notrack, Int_t paddleid, Int_t adcfl, Float_t &PadEdx, Int_t &SatWarning){
540    
541    //  ToFTrkVar *trk = GetToFTrkVar(notrack);
542        ToFTrkVar *trk = GetToFStoredTrack(notrack); //Elena 2015
543      this->GetdEdxPaddle(trk, paddleid, adcfl, PadEdx, SatWarning);
544      
545    };
546    
547    //
548    //  wm Nov 08 revision - saturation values included
549    /// gf Apr 07
550    /**
551     * Method to get the dEdx from a given ToF paddle.
552     * If two PMTs are good, the mean dEdx of both PMTs is taken, otherwise
553     * just the dEdx of the "good" PMT. If both PMTs are above saturation => dEdx=1000
554     * @param notrack Track Number
555     * @param Paddle index (0,1,...,23).
556     * @param adcflag in the paddle (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )
557     * @param PadEdx dEdx from a given ToF paddle
558     * @param SatWarning 1 if the PMT ios near saturation region (adcraw ~3000)
559     */
560    void ToFLevel2::GetdEdxPaddle(ToFTrkVar *trk, Int_t paddleid, Int_t adcfl, Float_t &PadEdx, Int_t &SatWarning){
561    
562      /*
563        Float_t  PMTsat[48] = {
564        3162.14, 3165.48, 3153.85, 3085.73, 3089.65, 3107.64, 3097.52, 3078.37,
565        3130.05, 3087.07, 3112.22, 3102.92, 3080.58, 3092.55, 3087.94, 3125.03,
566        3094.09, 3143.16, 3125.51, 3181.27, 3092.09, 3124.98, 3069.3, 3095.53,
567        3097.11, 3133.53, 3114.73, 3113.01, 3091.19, 3097.99, 3033.84, 3134.98,
568        3081.37, 3111.04, 3066.77, 3108.17, 3133, 3111.06, 3052.52, 3140.66,
569        3106.33, 3094.85, 3150.85, 3118.8, 3096.24, 3118.47,3111.36, 3117.11 } ;
570      */
571    
572      // new values from Napoli dec 2008
573      Float_t  PMTsat[48] = {
574        3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
575        3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
576        3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
577        3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
578        3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
579        3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };
580    
581      for (Int_t i=0; i<48;i++) PMTsat[i] = PMTsat[i] - 5.;  // safety margin
582    
583    
584      PadEdx = 0.;
585      //  SatWarning = 1000;
586      SatWarning = 0;   // 0=good, increase for each bad PMT
587    
588      Float_t dEdx[48] = {0};
589      Int_t pmt_id = -1;
590      Float_t adcraw[48];
591      //
592      if(!trk)cout << "ToFLevel2::GetdEdxPaddle(...) ---> NULL ToFTrkVar obj "<<endl;
593      if(!trk) return; //ELENA
594      //
595    
596      Int_t pmtleft=-1;
597      Int_t pmtright=-1;
598      GetPaddlePMT(paddleid, pmtleft, pmtright);
599    
600      adcraw[pmtleft] = 4095;
601      adcraw[pmtright] = 4095;
602    
603      
604      for (Int_t jj=0; jj<npmt(); jj++){
605        
606        ToFPMT *pmt = GetToFPMT(jj);
607        if(!pmt)break; //ELENA
608        
609        pmt_id = pmt->pmt_id;
610        if(pmt_id==pmtleft){
611          adcraw[pmtleft] = pmt->adc;
612        }
613        
614        if(pmt_id==pmtright){
615          adcraw[pmtright] = pmt->adc;
616        }
617      }
618    
619      
620      for (Int_t i=0; i<trk->npmtadc; i++){
621    
622        if((trk->adcflag).At(i)==0 || adcfl==100){
623          if((trk->pmtadc).At(i) == pmtleft)dEdx[pmtleft] = (trk->dedx).At(i);
624          if((trk->pmtadc).At(i) == pmtright)dEdx[pmtright] = (trk->dedx).At(i);
625        }else{
626          if((trk->pmtadc).At(i) == pmtleft)dEdx[pmtleft] = 0.;
627          if((trk->pmtadc).At(i) == pmtright)dEdx[pmtright] = 0.;
628        }
629      }
630    
631    
632      //  if( adcraw[pmtleft] >3000 || adcraw[pmtright] >3000)SatWarning=1;  //old version
633    
634      // Increase SatWarning Counter for each PMT>Sat
635      if( adcraw[pmtleft] > PMTsat[pmtleft])SatWarning++;  
636      if( adcraw[pmtright] > PMTsat[pmtright])SatWarning++;
637    
638      // if ADC  > sat set dEdx=1000
639      if( adcraw[pmtleft] > PMTsat[pmtleft]) dEdx[pmtleft] = 1000.;
640      if( adcraw[pmtright] > PMTsat[pmtright]) dEdx[pmtright] = 1000. ;
641    
642      // if two PMT are good, take mean dEdx, otherwise only the good dEdx
643      if(dEdx[pmtleft]<1000 && dEdx[pmtright]<1000) PadEdx = (dEdx[pmtleft]+dEdx[pmtright])*0.5;
644      if(dEdx[pmtleft]==1000 && dEdx[pmtright]<1000) PadEdx = dEdx[pmtright];  
645      if(dEdx[pmtleft]<1000 && dEdx[pmtright]==1000) PadEdx = dEdx[pmtleft];
646      
647    };
648    
649    // gf Apr 07
650    
651    /**
652     * Method to get the PMT name (like "S11_1A") if the PMT_ID is given.
653     * Indexes of corresponding  plane, paddle and  pmt are also given as output.
654     * @param ind  PMT_ID (0 - 47)
655     * @param iplane plane index (0 - 5)
656     * @param ipaddle paddle index (relative to the plane)
657     * @param ipmt pmt index (0(A), 1(B))
658     */
659    TString ToFLevel2::GetPMTName(Int_t ind, Int_t &iplane, Int_t &ipaddle,Int_t &ipmt){
660      
661      TString pmtname = " ";
662      
663      static const TString photoS[48] = {
664        "S11_1A", "S11_1B", "S11_2A", "S11_2B", "S11_3A", "S11_3B", "S11_4A",
665        "S11_4B",
666        "S11_5A", "S11_5B", "S11_6A", "S11_6B", "S11_7A", "S11_7B", "S11_8A",
667        "S11_8B",
668        "S12_1A", "S12_1B", "S12_2A", "S12_2B", "S12_3A", "S12_3B", "S12_4A",
669        "S12_4B", "S12_5A",  "S12_5B", "S12_6A", "S12_6B",
670        "S21_1A", "S21_1B", "S21_2A", "S21_2B",
671        "S22_1A", "S22_1B", "S22_2A", "S22_2B",
672        "S31_1A", "S31_1B", "S31_2A", "S31_2B", "S31_3A", "S31_3B",
673        "S32_1A", "S32_1B", "S32_2A", "S32_2B", "S32_3A", "S32_3B"
674      };
675      
676      
677      pmtname = photoS[ind].Data();
678      
679      TString ss = pmtname(1,2);
680      iplane  = (int)(atoi(ss.Data())/10)*2-3+atoi(ss.Data())%10;
681      ss = pmtname(4);
682      ipaddle = atoi(ss.Data())-1 ;
683      if( pmtname.Contains("A") )ipmt=0;
684      if( pmtname.Contains("B") )ipmt=1;
685      
686      return pmtname;
687    };
688    /**
689     * Method to get the PMT name (like "S11_1A") if the PMT_ID is given
690     * @param ind  PMT_ID (0 - 47)
691     */
692    TString ToFLevel2::GetPMTName(Int_t ind){
693    
694      Int_t iplane  = -1;
695      Int_t ipaddle = -1;
696      Int_t ipmt    = -1;
697      return GetPMTName(ind,iplane,ipaddle,ipmt);
698      
699    };
700    
701    // wm jun 08
702    Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane){
703    return GetPaddleIdOfTrack(xtr ,ytr ,plane, 0.4);
704    }
705    
706    // gf Apr 07
707    Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane, Float_t margin){
708      
709      Double_t xt,yt,xl,xh,yl,yh;
710      
711      Float_t tof11_x[8] = {-17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85};
712      Float_t tof12_y[6] = { -13.75,-8.25,-2.75,2.75,8.25,13.75};
713      Float_t tof21_y[2] = { 3.75,-3.75};
714      Float_t tof22_x[2] = { -4.5,4.5};
715      Float_t tof31_x[3] = { -6.0,0.,6.0};
716      Float_t tof32_y[3] = { -5.0,0.0,5.0};
717      
718      //  S11 8 paddles  33.0 x 5.1 cm
719      //  S12 6 paddles  40.8 x 5.5 cm
720      //  S21 2 paddles  18.0 x 7.5 cm
721      //  S22 2 paddles  15.0 x 9.0 cm
722      //  S31 3 paddles  15.0 x 6.0 cm
723      //  S32 3 paddles  18.0 x 5.0 cm
724      
725      Int_t paddleidoftrack=-1;
726      //
727      
728      //--- S11 ------
729      
730      if(plane==0){
731        xt = xtr;
732        yt = ytr;
733        paddleidoftrack=-1;
734        yl = -33.0/2. ;
735        yh =  33.0/2. ;
736        if ((yt>yl)&&(yt<yh)) {
737          for (Int_t i1=0; i1<8;i1++){
738            xl = tof11_x[i1] - (5.1-margin)/2. ;
739            xh = tof11_x[i1] + (5.1-margin)/2. ;
740            if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
741          }
742        }
743      }
744      //      cout<<"S11  "<<paddleidoftrack[0]<<"\n";
745      
746      //--- S12 -------
747      if(plane==1){
748        xt = xtr;
749        yt = ytr;
750        paddleidoftrack=-1;
751        xl = -40.8/2. ;
752        xh =  40.8/2. ;
753        
754        if ((xt>xl)&&(xt<xh)) {
755          for (Int_t i1=0; i1<6;i1++){
756            yl = tof12_y[i1] - (5.5-margin)/2. ;
757            yh = tof12_y[i1] + (5.5-margin)/2. ;
758            if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;
759          }
760        }
761      }
762      
763      //--- S21 ------
764    
765      if(plane==2){
766        xt = xtr;
767        yt = ytr;
768        paddleidoftrack=-1;
769        xl = -18./2. ;
770        xh =  18./2. ;
771        
772        if ((xt>xl)&&(xt<xh)) {
773          for (Int_t i1=0; i1<2;i1++){
774            yl = tof21_y[i1] - (7.5-margin)/2. ;
775            yh = tof21_y[i1] + (7.5-margin)/2. ;
776            if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;
777          }
778        }
779      }
780      
781      //--- S22 ------
782      if(plane==3){
783        xt = xtr;
784        yt = ytr;
785        paddleidoftrack=-1;
786        yl = -15./2. ;
787        yh =  15./2. ;
788        
789        if ((yt>yl)&&(yt<yh)) {
790          for (Int_t i1=0; i1<2;i1++){
791            xl = tof22_x[i1] - (9.0-margin)/2. ;
792            xh = tof22_x[i1] + (9.0-margin)/2. ;
793            if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
794          }
795        }
796      }  
797      
798      //--- S31 ------
799      if(plane==4){
800        xt = xtr;
801        yt = ytr;
802        paddleidoftrack=-1;
803        yl = -15.0/2. ;
804        yh =  15.0/2. ;
805        
806        if ((yt>yl)&&(yt<yh)) {
807          for (Int_t i1=0; i1<3;i1++){
808            xl = tof31_x[i1] - (6.0-margin)/2. ;
809            xh = tof31_x[i1] + (6.0-margin)/2. ;
810            if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
811          }
812        }
813      }  
814      
815      //---  S32 ------
816      if(plane==5){
817        xt = xtr;
818        yt = ytr;
819        paddleidoftrack=-1;
820        xl = -18.0/2. ;
821        xh =  18.0/2. ;
822        
823        if ((xt>xl)&&(xt<xh)) {
824          for (Int_t i1=0; i1<3;i1++){
825            yl = tof32_y[i1] - (5.0-margin)/2. ;
826            yh = tof32_y[i1] + (5.0-margin)/2. ;
827            if ((yt>yl)&&(yt<yh)) paddleidoftrack=i1;
828          }
829        }
830      }
831      
832      return paddleidoftrack;
833    
834    }  
835    
836    //
837    
838    // gf Apr 07
839    
840    void ToFLevel2::GetPMTPaddle(Int_t pmt_id, Int_t &plane, Int_t &paddle){
841      
842      plane = GetPlaneIndex(pmt_id);
843    
844      if(plane == 0){
845        if(pmt_id==0 || pmt_id==1)paddle=0;
846        if(pmt_id==2 || pmt_id==3)paddle=1;
847        if(pmt_id==4 || pmt_id==5)paddle=2;
848        if(pmt_id==6 || pmt_id==7)paddle=3;
849        if(pmt_id==8 || pmt_id==9)paddle=4;
850        if(pmt_id==10 || pmt_id==11)paddle=5;
851        if(pmt_id==12 || pmt_id==13)paddle=6;
852        if(pmt_id==14 || pmt_id==15)paddle=7;
853      }
854      
855      if(plane == 1){
856        if(pmt_id==16 || pmt_id==17)paddle=0;
857        if(pmt_id==18 || pmt_id==19)paddle=1;
858        if(pmt_id==20 || pmt_id==21)paddle=2;
859        if(pmt_id==22 || pmt_id==23)paddle=3;
860        if(pmt_id==24 || pmt_id==25)paddle=4;
861        if(pmt_id==26 || pmt_id==27)paddle=5;
862      }
863      
864      if(plane == 2){
865        if(pmt_id==28 || pmt_id==29)paddle=0;
866        if(pmt_id==30 || pmt_id==31)paddle=1;
867      }
868      
869      if(plane == 3){
870        if(pmt_id==32 || pmt_id==33)paddle=0;
871        if(pmt_id==34 || pmt_id==35)paddle=1;
872      }
873      
874      if(plane == 4){
875        if(pmt_id==36 || pmt_id==37)paddle=0;
876        if(pmt_id==38 || pmt_id==39)paddle=1;
877        if(pmt_id==40 || pmt_id==41)paddle=2;
878      }
879      
880      if(plane == 5){
881        if(pmt_id==42 || pmt_id==43)paddle=0;
882        if(pmt_id==44 || pmt_id==45)paddle=1;
883        if(pmt_id==46 || pmt_id==47)paddle=2;
884      }
885      return;
886    }
887    
888    //
889    
890    // gf Apr 07
891    
892    void ToFLevel2::GetPaddlePMT(Int_t paddle, Int_t &pmtleft, Int_t &pmtright){
893      pmtleft=paddle*2;
894      pmtright= pmtleft+1;  
895      return;
896    }
897    
898    //
899    
900    
901    
902    // // gf Apr 07
903    
904    void ToFLevel2::GetPaddleGeometry(Int_t plane, Int_t paddle, Float_t &xleft, Float_t &xright, Float_t &yleft, Float_t &yright){
905      
906      Int_t i1;
907    
908      Float_t tof11_x[8] = {-17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85};
909      Float_t tof12_y[6] = { -13.75,-8.25,-2.75,2.75,8.25,13.75};
910      Float_t tof21_y[2] = { 3.75,-3.75};
911      Float_t tof22_x[2] = { -4.5,4.5};
912      Float_t tof31_x[3] = { -6.0,0.,6.0};
913      Float_t tof32_y[3] = { -5.0,0.0,5.0};
914            
915      //  S11 8 paddles  33.0 x 5.1 cm
916      //  S12 6 paddles  40.8 x 5.5 cm
917      //  S21 2 paddles  18.0 x 7.5 cm
918      //  S22 2 paddles  15.0 x 9.0 cm
919      //  S31 3 paddles  15.0 x 6.0 cm
920      //  S32 3 paddles  18.0 x 5.0 cm
921    
922      if(plane==0)
923        {
924          for (i1=0; i1<8;i1++){
925            if(i1 == paddle){
926              xleft = tof11_x[i1] - 5.1/2.;
927              xright = tof11_x[i1] + 5.1/2.;
928              yleft = -33.0/2.;
929              yright = 33.0/2.;
930            }
931          }
932        }
933      
934      if(plane==1)
935        {
936          for (i1=0; i1<6;i1++){
937            if(i1 == paddle){
938              xleft = -40.8/2.;
939              xright = 40.8/2.;
940              yleft = tof12_y[i1] - 5.5/2.;
941              yright = tof12_y[i1] + 5.5/2.;
942            }
943          }
944        }
945    
946      if(plane==2)
947        {
948          for (i1=0; i1<2;i1++){
949            if(i1 == paddle){
950              xleft =  -18./2.;
951              xright = 18./2.;
952              yleft = tof21_y[i1] - 7.5/2.;
953              yright = tof21_y[i1] + 7.5/2.;
954            }
955          }
956        }
957      
958      if(plane==3)
959        {
960          for (i1=0; i1<2;i1++){
961            if(i1 == paddle){
962              xleft = tof22_x[i1] - 9.0/2.;
963              xright = tof22_x[i1] + 9.0/2.;
964              yleft = -15./2.;
965              yright = 15./2.;
966            }
967          }
968        }
969    
970    
971      if(plane==4)
972        {
973          for (i1=0; i1<3;i1++){
974            if(i1 == paddle){
975              xleft = tof31_x[i1] - 6.0/2.;
976              xright = tof31_x[i1] + 6.0/2.;
977              yleft = -15./2.;
978              yright = 15./2.;
979            }
980          }
981        }
982    
983      if(plane==5)
984        {
985          for (i1=0; i1<3;i1++){
986            if(i1 == paddle){
987              xleft = -18.0/2.;
988              xright = 18.0/2.;
989              yleft = tof32_y[i1] - 5.0/2.;
990              yright = tof32_y[i1] + 5.0/2.;
991            }
992          }
993        }
994      return;
995    }
996    
997    // gf Apr 07
998    /**
999     * Method to get the paddle index (0,...23) if the plane ID and the paddle id in the plane is given.
1000     * This method is the
1001     * "reverse" of method "GetPaddlePlane"
1002     * @param plane    (0 - 5)
1003     * @param paddle   (plane=0, paddle = 0,...5)
1004     * @param padid    (0 - 23)
1005     */
1006    Int_t ToFLevel2::GetPaddleid(Int_t plane, Int_t paddle)
1007    {
1008      Int_t padid=-1;
1009      Int_t pads[6]={8,6,2,2,3,3};
1010    
1011      int somma=0;
1012      int np=plane;
1013      for(Int_t j=0; j<np; j++){
1014        somma+=pads[j];
1015      }
1016      padid=paddle+somma;
1017      return padid;
1018    
1019    }
1020    
1021    
1022    // gf Apr 07
1023    /**
1024     * Method to get the plane ID and the paddle id in the plane if the paddle index (0,...23) is given.
1025     * This method is the
1026     * "reverse" of method "GetPaddleid"
1027     * @param pad      (0 - 23)
1028     * @param plane    (0 - 5)
1029     * @param paddle   (plane=0, paddle = 0,...5)
1030     */
1031    void ToFLevel2::GetPaddlePlane(Int_t pad, Int_t &plane, Int_t &paddle)
1032    {
1033    
1034      Int_t pads11=8;
1035      Int_t pads12=6;
1036      Int_t pads21=2;
1037      Int_t pads22=2;
1038      Int_t pads31=3;
1039      // Int_t pads32=3;
1040    
1041      if(pad<8){
1042        plane=0;
1043        paddle=pad;
1044        return;
1045      }
1046    
1047      if((7<pad)&&(pad<14)){
1048        plane=1;
1049        paddle=pad-pads11;
1050        return;
1051      }
1052      
1053      if((13<pad)&&(pad<16)){
1054        plane=2;
1055        paddle=pad-pads11-pads12;
1056        return;
1057      }
1058    
1059      if((15<pad)&&(pad<18)){
1060        plane=3;
1061        paddle=pad-pads11-pads12-pads21;
1062        return;
1063      }
1064    
1065      if((17<pad)&&(pad<21)){
1066        plane=4;
1067        paddle=pad-pads11-pads12-pads21-pads22;
1068        return;
1069      }
1070    
1071      if((20<pad)&&(pad<24)){
1072        plane=5;
1073        paddle=pad-pads11-pads12-pads21-pads22-pads31;
1074        return;
1075      }  
1076    
1077    }
1078    
1079    
1080    Int_t ToFLevel2::GetNPaddle(Int_t plane){
1081    
1082      Int_t npaddle=-1;
1083    
1084      Int_t pads11=8;
1085      Int_t pads12=6;
1086      Int_t pads21=2;
1087      Int_t pads22=2;
1088      Int_t pads31=3;
1089      Int_t pads32=3;
1090    
1091      if(plane==0)npaddle=pads11;
1092      if(plane==1)npaddle=pads12;
1093      if(plane==2)npaddle=pads21;
1094      if(plane==3)npaddle=pads22;
1095      if(plane==4)npaddle=pads31;
1096      if(plane==5)npaddle=pads32;
1097    
1098      return npaddle;
1099    
1100    }
1101    
1102    
1103    
1104    /// wm feb 08
1105    
1106    /**
1107     * Method to calculate Beta from the 12 single measurements
1108     * we check the individual weights for artificial TDC values, then calculate
1109     * am mean beta for the first time. In a second step we loop again through
1110     * the single measurements, checking for the residual from the mean
1111     * The cut on the residual reject measurements > "x"-sigma. A chi2 value is
1112     * calculated, furthermore a "quality" value by adding the weights which
1113     * are finally used. If all measurements are taken, "quality" will be = 22.47.
1114     * A chi2 cut around 3-4 and a quality-cut > 20 is needed for clean beta
1115     * measurements like antiprotons etc.
1116     * The Level2 output is derived in the fortran routines using: 10.,10.,20.
1117     * @param notrack Track Number
1118     * @param cut on residual: difference between single measurement and mean
1119     * @param cut on "quality"
1120     * @param cut on chi2
1121     */
1122    
1123    
1124    Float_t ToFTrkVar::CalcBeta( Float_t resmax, Float_t qualitycut, Float_t chi2cut){
1125    
1126    
1127      Float_t bxx = 100.;
1128      //
1129      ToFTrkVar *trk = this;
1130    
1131    
1132      Float_t chi2,xhelp,beta_mean;
1133      Float_t w_i[12],quality,sw,sxw,res,betachi,beta_mean_inv;
1134      Float_t b[12],tdcfl;
1135      Int_t  pmt_id,pmt_plane;
1136    
1137      for (Int_t i=0; i<12; i++){
1138        b[i] = trk->beta[i];
1139                                  }
1140          
1141    
1142    //========================================================================
1143    //---  Find out ToF layers with artificial TDC values & fill vector    ---
1144    //========================================================================
1145    
1146    Float_t  w_il[6];
1147    
1148         for (Int_t jj=0; jj<6;jj++) {
1149             w_il[jj] = 1000.;
1150                                     }
1151    
1152    
1153      for (Int_t i=0; i<trk->npmttdc; i++){
1154        //
1155        pmt_id = (trk->pmttdc).At(i);
1156        pmt_plane = ToFLevel2::GetPlaneIndex(pmt_id);
1157        tdcfl = (trk->tdcflag).At(i);
1158        if (w_il[pmt_plane] != 1.) w_il[pmt_plane] = tdcfl; //tdcflag
1159                                         };
1160      
1161    //========================================================================
1162    //---  Set weights for the 12 measurements using information for top and bottom:
1163    //---  if no measurements: weight = set to very high value=> not used
1164    //---  top or bottom artificial: weight*sqrt(2)
1165    //---  top and bottom artificial: weight*sqrt(2)*sqrt(2)
1166    //========================================================================
1167    
1168    Int_t itop[12] = {0,0,1,1,2,2,3,3,0,0,1,1};
1169    Int_t ibot[12] = {4,5,4,5,4,5,4,5,2,3,2,3};
1170    
1171         xhelp= 1E09;
1172      
1173         for (Int_t jj=0; jj<12;jj++) {
1174         if (jj<4)           xhelp = 0.11;    // S1-S3
1175         if ((jj>3)&&(jj<8)) xhelp = 0.18;    // S2-S3
1176         if (jj>7)           xhelp = 0.28;    // S1-S2
1177         if ((w_il[itop[jj]] == 1000.) && (w_il[ibot[jj]] == 1000.)) xhelp = 1E09;
1178         if ((w_il[itop[jj]] == 1) || (w_il[ibot[jj]] == 1.)) xhelp = xhelp*1.414 ;
1179         if ((w_il[itop[jj]] == 1) && (w_il[ibot[jj]] == 1.)) xhelp = xhelp*2. ;
1180    
1181         w_i[jj] = 1./xhelp;
1182                                      }
1183    
1184    
1185    //========================================================================
1186    //--- Calculate mean beta for the first time -----------------------------
1187    //--- We are using "1/beta" since its error is gaussian ------------------
1188    //========================================================================
1189    
1190          Int_t icount=0;
1191          sw=0.;
1192          sxw=0.;
1193          beta_mean=100.;
1194    
1195              for (Int_t jj=0; jj<12;jj++){
1196            if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.))
1197             {
1198                icount= icount+1;
1199                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1200                sw =sw + w_i[jj]*w_i[jj] ;
1201    
1202             }
1203             }
1204    
1205          if (icount>0) beta_mean=1./(sxw/sw);
1206          beta_mean_inv = 1./beta_mean;
1207    
1208    //========================================================================
1209    //--- Calculate beta for the second time, use residuals of the single
1210    //--- measurements to get a chi2 value
1211    //========================================================================
1212    
1213          icount=0;
1214          sw=0.;
1215          sxw=0.;
1216          betachi = 100.;
1217          chi2 = 0.;
1218          quality=0.;
1219    
1220    
1221              for (Int_t jj=0; jj<12;jj++){
1222           if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.)&&(w_i[jj]>0.01)) {
1223                res = beta_mean_inv - (1./b[jj]) ;
1224                if (fabs(res*w_i[jj])<resmax)          {;
1225                chi2 = chi2 + pow((res*w_i[jj]),2) ;
1226                icount= icount+1;
1227                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1228                sw =sw + w_i[jj]*w_i[jj] ;
1229                                                   }
1230                                                                            }
1231                                          }
1232          quality = sqrt(sw) ;
1233    
1234          if (icount==0) chi2 = 1000.;
1235          if (icount>0) chi2 = chi2/(icount) ;
1236          if (icount>0) betachi=1./(sxw/sw);
1237    
1238       bxx = 100.;
1239       if ((chi2 < chi2cut)&&(quality>qualitycut)) bxx = betachi;
1240      //
1241      return(bxx);
1242    };
1243    ////////////////////////////////////////////////////
1244    ////////////////////////////////////////////////////
1245    /**
1246     * See ToFTrkVar::CalcBeta(Float_t,Float_t, Float_t).
1247     */
1248    Float_t ToFLevel2::CalcBeta(Int_t notrack, Float_t resmax, Float_t qualitycut, Float_t chi2cut){
1249    
1250    //  cout<<" in CalcBeta "<<resmax<<" "<<chi2cut<<" "<<qualitycut<<endl;
1251    
1252    //  ToFTrkVar *trk = GetToFTrkVar(notrack); //wrong!
1253        ToFTrkVar *trk = GetToFStoredTrack(notrack);//Elena Apr 2015
1254      if(!trk) return 0; //ELENA
1255    
1256      return trk->CalcBeta(resmax,qualitycut,chi2cut);
1257    
1258    };
1259    
1260    
1261    ////////////////////////////////////////////////////
1262    ////////////////////////////////////////////////////
1263    
1264    
1265    /**
1266     * Fills a struct cToFLevel2 with values from a ToFLevel2 object (to put data into a F77 common).
1267     */
1268    void ToFLevel2::GetLevel2Struct(cToFLevel2 *l2) const{
1269    
1270      for(Int_t i=0;i<6;i++)
1271        l2->tof_j_flag[i]=tof_j_flag[i];
1272    
1273      if(ToFTrk){ //ELENA
1274          l2->ntoftrk = ToFTrk->GetEntries();
1275          for(Int_t j=0;j<l2->ntoftrk;j++){
1276              l2->toftrkseqno[j]= ((ToFTrkVar*)ToFTrk->At(j))->trkseqno;
1277              l2->npmttdc[j]= ((ToFTrkVar*)ToFTrk->At(j))->npmttdc;
1278              for(Int_t i=0;i<l2->npmttdc[j];i++){
1279                  l2->pmttdc[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->pmttdc.At(i);
1280                  l2->tdcflag[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->tdcflag.At(i); // gf: 30 Nov 2006
1281              }
1282              for(Int_t i=0;i<13;i++)
1283                  l2->beta[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->beta[i];
1284              
1285              l2->npmtadc[j]= ((ToFTrkVar*)ToFTrk->At(j))->npmtadc;
1286              for(Int_t i=0;i<l2->npmtadc[j];i++){
1287                  l2->pmtadc[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->pmtadc.At(i);
1288                  l2->adcflag[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->adcflag.At(i); // gf: 30 Nov 2006
1289                  l2->dedx[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->dedx.At(i);
1290              }
1291              for(Int_t i=0;i<3;i++){
1292                  l2->xtofpos[i][j]=((ToFTrkVar*)ToFTrk->At(j))->xtofpos[i];
1293                  l2->ytofpos[i][j]=((ToFTrkVar*)ToFTrk->At(j))->ytofpos[i];
1294              }
1295              for(Int_t i=0;i<6;i++){
1296                  l2->xtr_tof[i][j]=((ToFTrkVar*)ToFTrk->At(j))->xtr_tof[i];
1297                  l2->ytr_tof[i][j]=((ToFTrkVar*)ToFTrk->At(j))->ytr_tof[i];
1298              }
1299          }
1300      } //ELENA
1301        
1302      if(PMT){ //ELENA
1303          l2->npmt = PMT->GetEntries();
1304          for(Int_t j=0;j<l2->npmt;j++){
1305              l2->pmt_id[j] = ((ToFPMT*)PMT->At(j))->pmt_id;
1306              l2->adc[j] =((ToFPMT*)PMT->At(j))->adc;
1307              l2->tdc_tw[j] =((ToFPMT*)PMT->At(j))->tdc_tw;
1308          }
1309      } //ELENA
1310    }
1311    
1312    
1313    //
1314    // Reprocessing tool // Emiliano 08/04/07
1315    //
1316    Int_t ToFLevel2::Process(TrkLevel2 *trk, TrigLevel2 *trg, GL_RUN *run, OrbitalInfo *orb, Bool_t force){
1317      //
1318      // Copiare qui qualcosa di simile a calonuclei per evitare di riprocessare sempre tutto
1319      //
1320      printf("\n\n\n ERROR: NOT IMPLEMENTED ANYMORE, write Emiliano if you need this method (Emiliano.Mocchiutti@ts.infn.it) \n\n\n");
1321      return(-1);
1322      //   //
1323      //   // structures to communicate with F77
1324      //   //
1325      //   extern struct ToFInput  tofinput_;
1326    //   extern struct ToFOutput tofoutput_;
1327    //   //
1328    //   // DB connection
1329    //   //
1330    //   TString host;
1331    //   TString user;
1332    //   TString psw;
1333    //   const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
1334    //   const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
1335    //   const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
1336    //   if ( !pamdbhost ) pamdbhost = "";
1337    //   if ( !pamdbuser ) pamdbuser = "";
1338    //   if ( !pamdbpsw ) pamdbpsw = "";
1339    //   if ( strcmp(pamdbhost,"") ) host = pamdbhost;
1340    //   if ( strcmp(pamdbuser,"") ) user = pamdbuser;
1341    //   if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
1342    //   //
1343    //   //
1344    //   TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
1345    //   if ( !dbc->IsConnected() ) return 1;
1346    //   stringstream myquery;
1347    //   myquery.str("");
1348    //   myquery << "SET time_zone='+0:00';";
1349    //   dbc->Query(myquery.str().c_str());
1350    //   delete dbc->Query("SET sql_mode = 'NO_UNSIGNED_SUBTRACTION';");
1351    //   GL_PARAM *glparam = new GL_PARAM();
1352    //   glparam->Query_GL_PARAM(1,1,dbc); // parameters stored in DB in GL_PRAM table
1353    //   trk->LoadField(glparam->PATH+glparam->NAME);
1354    //   //
1355    //   Bool_t defcal = true;
1356    //   Int_t error=glparam->Query_GL_PARAM(run->RUNHEADER_TIME,201,dbc); // parameters stored in DB in GL_PRAM table
1357    //   if ( error<0 ) {
1358    //     return(1);
1359    //   };
1360    //   printf(" Reading ToF parameter file: %s \n",(glparam->PATH+glparam->NAME).Data());
1361    //   if ( (UInt_t)glparam->TO_TIME != (UInt_t)4294967295UL ) defcal = false;
1362    //   //
1363    //   Int_t nlen = (Int_t)(glparam->PATH+glparam->NAME).Length();
1364    //   rdtofcal((char *)(glparam->PATH+glparam->NAME).Data(),&nlen);
1365    //   //
1366    //   Int_t adc[4][12];
1367    //   Int_t tdc[4][12];
1368    //   Float_t tdcc[4][12];
1369    //   //
1370    //   // process tof data
1371    //   //
1372    //   for (Int_t hh=0; hh<12;hh++){
1373    //     for (Int_t kk=0; kk<4;kk++){
1374    //            adc[kk][hh] = 4095;
1375    //            tdc[kk][hh] = 4095;
1376    //            tdcc[kk][hh] = 4095.;
1377    //            tofinput_.adc[hh][kk] = 4095;
1378    //            tofinput_.tdc[hh][kk] = 4095;
1379    //     };
1380    //   };
1381    //   Int_t ntrkentry = 0;
1382    //   Int_t npmtentry = 0;
1383    //   Int_t gg = 0;
1384    //   Int_t hh = 0;
1385    //   Int_t adcf[48];
1386    //   memset(adcf, 0, 48*sizeof(Int_t));
1387    //   Int_t tdcf[48];
1388    //   memset(tdcf, 0, 48*sizeof(Int_t));
1389    //   for (Int_t pm=0; pm < this->ntrk() ; pm++){
1390    //      ToFTrkVar *ttf = this->GetToFTrkVar(pm);
1391    //      for ( Int_t nc=0; nc < ttf->npmttdc; nc++){
1392    //             if ( (ttf->tdcflag).At(nc) != 0 ) tdcf[(ttf->pmttdc).At(nc)] = 1;
1393    //      };
1394    //      for ( Int_t nc=0; nc < ttf->npmtadc; nc++){
1395    //             if ( (ttf->adcflag).At(nc) != 0 ) adcf[(ttf->pmtadc).At(nc)] = 1;
1396    //      };
1397    //   };
1398    //   //
1399    //   for (Int_t pm=0; pm < this->npmt() ; pm++){
1400    //      ToFPMT *pmt = this->GetToFPMT(pm);
1401    //      this->GetPMTIndex(pmt->pmt_id, gg, hh);
1402    //      if ( adcf[pmt->pmt_id] == 0 ){
1403    //              tofinput_.adc[gg][hh] = (int)pmt->adc;
1404    //              adc[hh][gg] = (int)pmt->adc;
1405    //      };
1406    //      if ( tdcf[pmt->pmt_id] == 0 ){
1407    //              tofinput_.tdc[gg][hh] = (int)pmt->tdc;
1408    //              tdc[hh][gg] = (int)pmt->tdc;
1409    //      };
1410    //      tdcc[hh][gg] = (float)pmt->tdc_tw;
1411    //      // Int_t pppid = this->GetPMTid(hh,gg);
1412    //      //      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);
1413    //   };
1414    //   //
1415    //   Int_t unpackError = this->unpackError;
1416    //   //
1417    //   for (Int_t hh=0; hh<5;hh++){
1418    //      tofinput_.patterntrig[hh]=trg->patterntrig[hh];
1419    //   };
1420    //   //
1421    //   this->Clear();
1422    //   //
1423    //       Int_t pmt_id = 0;
1424    //       ToFPMT *t_pmt = new ToFPMT();
1425    //       if(!(this->PMT)) this->PMT = new TClonesArray("ToFPMT",12); //ELENA
1426    //       TClonesArray &tpmt = *this->PMT;
1427    //       ToFTrkVar *t_tof = new ToFTrkVar();
1428    //       if(!(this->ToFTrk)) this->ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
1429    //       TClonesArray &t = *this->ToFTrk;
1430    //       //
1431    //       //
1432    //       // Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract common and track-related  variables.
1433    //       //
1434    //       npmtentry = 0;
1435    //       //
1436    //       ntrkentry = 0;
1437    //       //
1438    //       // Calculate tracks informations from ToF alone
1439    //       //
1440    //       tofl2com();
1441    //       //
1442    //       memcpy(this->tof_j_flag,tofoutput_.tof_j_flag,6*sizeof(Int_t));
1443    //       //
1444    //       t_tof->trkseqno = -1;
1445    //       //
1446    //       // and now we must copy from the output structure to the level2 class:
1447    //       //
1448    //       t_tof->npmttdc = 0;
1449    //       //
1450    //       for (Int_t hh=0; hh<12;hh++){
1451    //         for (Int_t kk=0; kk<4;kk++){
1452    //           if ( tofoutput_.tofmask[hh][kk] != 0 ){
1453    //             pmt_id = this->GetPMTid(kk,hh);
1454    //             t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
1455    //             t_tof->tdcflag.AddAt(tofoutput_.tdcflagtof[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
1456    //             t_tof->npmttdc++;
1457    //           };
1458    //         };
1459    //       };
1460    //       for (Int_t kk=0; kk<13;kk++){
1461    //         t_tof->beta[kk] = tofoutput_.betatof_a[kk];
1462    //       }
1463    //       //
1464    //       t_tof->npmtadc = 0;
1465    //       for (Int_t hh=0; hh<12;hh++){
1466    //         for (Int_t kk=0; kk<4;kk++){
1467    //           if ( tofoutput_.adctof_c[hh][kk] < 1000 ){
1468    //             t_tof->dedx.AddAt(tofoutput_.adctof_c[hh][kk],t_tof->npmtadc);
1469    //             pmt_id = this->GetPMTid(kk,hh);
1470    //             t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
1471    //             t_tof->adcflag.AddAt(tofoutput_.adcflagtof[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
1472    //             t_tof->npmtadc++;
1473    //           };
1474    //         };
1475    //       };
1476    //       //
1477    //       memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
1478    //       memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
1479    //       memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
1480    //       memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
1481    //       //
1482    //       new(t[ntrkentry]) ToFTrkVar(*t_tof);
1483    //       ntrkentry++;
1484    //       t_tof->Clear();
1485    //       //
1486    //       //
1487    //       //
1488    //       t_pmt->Clear();
1489    //       //
1490    //       for (Int_t hh=0; hh<12;hh++){
1491    //         for (Int_t kk=0; kk<4;kk++){
1492    //          // new WM
1493    //           if ( tofoutput_.tdc_c[hh][kk] < 4095 || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){
1494    // //          if ( tdcc[kk][hh] < 4095. || adc[kk][hh] < 4095  || tdc[kk][hh] < 4095 ){
1495    //             //
1496    //             t_pmt->pmt_id = this->GetPMTid(kk,hh);
1497    //             t_pmt->tdc_tw = tofoutput_.tdc_c[hh][kk];
1498    //             t_pmt->adc = (Float_t)adc[kk][hh];
1499    //             t_pmt->tdc = (Float_t)tdc[kk][hh];
1500    //             //
1501    //             new(tpmt[npmtentry]) ToFPMT(*t_pmt);
1502    //             npmtentry++;
1503    //             t_pmt->Clear();
1504    //           };
1505    //         };
1506    //       };
1507    //       //
1508    //       // Calculate track-related variables
1509    //       //
1510    //       if ( trk->ntrk() > 0 ){
1511    //         //
1512    //         // We have at least one track
1513    //         //
1514    //         //
1515    //         // Run over tracks
1516    //         //
1517    //         for(Int_t nt=0; nt < trk->ntrk(); nt++){
1518    //           //
1519    //           TrkTrack *ptt = trk->GetStoredTrack(nt);
1520    //           //
1521    //           // Copy the alpha vector in the input structure
1522    //           //
1523    //           for (Int_t e = 0; e < 5 ; e++){
1524    //             tofinput_.al_pp[e] = ptt->al[e];
1525    //           };
1526    //           //
1527    //           // Get tracker related variables for this track
1528    //           //
1529    //           toftrk();
1530    //           //
1531    //           // Copy values in the class from the structure (we need to use a temporary class to store variables).
1532    //           //
1533    //           t_tof->npmttdc = 0;
1534    //           for (Int_t hh=0; hh<12;hh++){
1535    //             for (Int_t kk=0; kk<4;kk++){
1536    //               if ( tofoutput_.tofmask[hh][kk] != 0 ){
1537    //                 pmt_id = this->GetPMTid(kk,hh);
1538    //                 t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc);
1539    //                 t_tof->tdcflag.AddAt(tofoutput_.tdcflag[hh][kk],t_tof->npmttdc); // gf: Jan 09/07
1540    //                 t_tof->npmttdc++;
1541    //               };
1542    //             };
1543    //           };
1544    //           for (Int_t kk=0; kk<13;kk++){
1545    //             t_tof->beta[kk] = tofoutput_.beta_a[kk];
1546    //           };
1547    //           //
1548    //           t_tof->npmtadc = 0;
1549    //           for (Int_t hh=0; hh<12;hh++){
1550    //             for (Int_t kk=0; kk<4;kk++){
1551    //               if ( tofoutput_.adc_c[hh][kk] < 1000 ){
1552    //                 t_tof->dedx.AddAt(tofoutput_.adc_c[hh][kk],t_tof->npmtadc);
1553    //                 pmt_id = this->GetPMTid(kk,hh);
1554    //                 t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc);
1555    //                 t_tof->adcflag.AddAt(tofoutput_.adcflag[hh][kk],t_tof->npmtadc); // gf: Jan 09/07
1556    //                 t_tof->npmtadc++;
1557    //               };
1558    //             };
1559    //           };
1560    //           //
1561    //           memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos));
1562    //           memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos));
1563    //           memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof));
1564    //           memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof));
1565    //           //
1566    //           // Store the tracker track number in order to be sure to have shyncronized data during analysis
1567    //           //
1568    //           t_tof->trkseqno = nt;
1569    //           //
1570    //           // create a new object for this event with track-related variables
1571    //           //
1572    //           new(t[ntrkentry]) ToFTrkVar(*t_tof);
1573    //           ntrkentry++;
1574    //           t_tof->Clear();
1575    //           //
1576    //         }; // loop on all the tracks
1577    //       //
1578    //       this->unpackError = unpackError;
1579    //       if ( defcal ){
1580    //         this->default_calib = 1;
1581    //       } else {
1582    //         this->default_calib = 0;
1583    //       };
1584    //};
1585    //  return(0);
1586    }
1587    
1588    bool ToFLevel2::bit(int decimal, char pos){
1589      return( (decimal>>pos)%2 );
1590    }
1591    
1592    bool ToFLevel2::checkPMT(TString givenpmt){
1593      TClonesArray* Pmt = this->PMT;
1594      //  printf(" ou %s entries %i \n",givenpmt.Data(),Pmt->GetEntries());
1595      for(int i=0; i<Pmt->GetEntries(); i++) {  
1596        ToFPMT* pmthit = (ToFPMT*)Pmt->At(i);
1597        TString pmtname = this->GetPMTName(pmthit->pmt_id);
1598        //    printf(" name %s \n",pmtname.Data());
1599        if ( !strcmp(pmtname.Data(),givenpmt.Data()) )
1600          return true;
1601      }
1602      //  printf(" PMT %s missing \n",givenpmt.Data());
1603      return false;
1604    }
1605    
1606    bool ToFLevel2::checkPMTpatternPMThit(TrigLevel2 *trg, int &pmtpattern, int &pmtnosignal){
1607      UInt_t *patterntrig = trg->patterntrig;
1608      pmtpattern = 0;
1609      pmtnosignal = 0;
1610      bool good = true;
1611      //S3
1612      if ( this->bit(patterntrig[2],0) ){ pmtpattern++;  if ( !this->checkPMT("S31_1A")){ pmtnosignal++; good = false;}}
1613      if ( this->bit(patterntrig[2],1) ){ pmtpattern++;  if ( !this->checkPMT("S31_2A")){ pmtnosignal++; good = false;}}
1614      if ( this->bit(patterntrig[2],2) ){ pmtpattern++;  if ( !this->checkPMT("S31_3A")){ pmtnosignal++; good = false;}}
1615      if ( this->bit(patterntrig[2],3) ){ pmtpattern++;  if ( !this->checkPMT("S31_1B")){ pmtnosignal++; good = false;}}
1616      if ( this->bit(patterntrig[2],4) ){ pmtpattern++;  if ( !this->checkPMT("S31_2B")){ pmtnosignal++; good = false;}}
1617      if ( this->bit(patterntrig[2],5) ){ pmtpattern++;  if ( !this->checkPMT("S31_3B")){ pmtnosignal++; good = false;}}      
1618      if ( this->bit(patterntrig[2],6) ){ pmtpattern++;  if ( !this->checkPMT("S32_1A")){ pmtnosignal++; good = false;}}
1619      if ( this->bit(patterntrig[2],7) ){ pmtpattern++;  if ( !this->checkPMT("S32_2A")){ pmtnosignal++; good = false;}}
1620      if ( this->bit(patterntrig[2],8) ){ pmtpattern++;  if ( !this->checkPMT("S32_3A")){ pmtnosignal++; good = false;}}
1621      if ( this->bit(patterntrig[2],9) ){ pmtpattern++;  if ( !this->checkPMT("S32_1B")){ pmtnosignal++; good = false;}}
1622      if ( this->bit(patterntrig[2],10) ){ pmtpattern++;  if ( !this->checkPMT("S32_2B")){ pmtnosignal++; good = false;}}
1623      if ( this->bit(patterntrig[2],11) ){ pmtpattern++;  if ( !this->checkPMT("S32_3B")){ pmtnosignal++; good = false;}}      
1624      //S2
1625      if ( this->bit(patterntrig[3],0) ){ pmtpattern++;  if ( !this->checkPMT("S21_1A")){ pmtnosignal++; good = false;}}
1626      if ( this->bit(patterntrig[3],1) ){ pmtpattern++;  if ( !this->checkPMT("S21_2A")){ pmtnosignal++; good = false;}}
1627      if ( this->bit(patterntrig[3],2) ){ pmtpattern++;  if ( !this->checkPMT("S21_1B")){ pmtnosignal++; good = false;}}
1628      if ( this->bit(patterntrig[3],3) ){ pmtpattern++;  if ( !this->checkPMT("S21_2B")){ pmtnosignal++; good = false;}}      
1629      if ( this->bit(patterntrig[3],4) ){ pmtpattern++;  if ( !this->checkPMT("S22_1A")){ pmtnosignal++; good = false;}}
1630      if ( this->bit(patterntrig[3],5) ){ pmtpattern++;  if ( !this->checkPMT("S22_2A")){ pmtnosignal++; good = false;}}
1631      if ( this->bit(patterntrig[3],6) ){ pmtpattern++;  if ( !this->checkPMT("S22_1B")){ pmtnosignal++; good = false;}}
1632      if ( this->bit(patterntrig[3],7) ){ pmtpattern++;  if ( !this->checkPMT("S22_2B")){ pmtnosignal++; good = false;}}      
1633      //S12
1634      if ( this->bit(patterntrig[4],0) ){ pmtpattern++;  if ( !this->checkPMT("S12_1A")){ pmtnosignal++; good = false;}}
1635      if ( this->bit(patterntrig[4],1) ){ pmtpattern++;  if ( !this->checkPMT("S12_2A")){ pmtnosignal++; good = false;}}
1636      if ( this->bit(patterntrig[4],2) ){ pmtpattern++;  if ( !this->checkPMT("S12_3A")){ pmtnosignal++; good = false;}}
1637      if ( this->bit(patterntrig[4],3) ){ pmtpattern++;  if ( !this->checkPMT("S12_4A")){ pmtnosignal++; good = false;}}
1638      if ( this->bit(patterntrig[4],4) ){ pmtpattern++;  if ( !this->checkPMT("S12_5A")){ pmtnosignal++; good = false;}}
1639      if ( this->bit(patterntrig[4],5) ){ pmtpattern++;  if ( !this->checkPMT("S12_6A")){ pmtnosignal++; good = false;}}      
1640      if ( this->bit(patterntrig[4],6) ){ pmtpattern++;  if ( !this->checkPMT("S12_1A")){ pmtnosignal++; good = false;}}
1641      if ( this->bit(patterntrig[4],7) ){ pmtpattern++;  if ( !this->checkPMT("S12_2A")){ pmtnosignal++; good = false;}}
1642      if ( this->bit(patterntrig[4],8) ){ pmtpattern++;  if ( !this->checkPMT("S12_3A")){ pmtnosignal++; good = false;}}
1643      if ( this->bit(patterntrig[4],9) ){ pmtpattern++;  if ( !this->checkPMT("S12_4B")){ pmtnosignal++; good = false;}}
1644      if ( this->bit(patterntrig[4],10) ){ pmtpattern++; if ( !this->checkPMT("S12_5B")){ pmtnosignal++; good = false;}}
1645      if ( this->bit(patterntrig[4],11) ){ pmtpattern++; if ( !this->checkPMT("S12_6B")){ pmtnosignal++; good = false;}}      
1646      //S11
1647      if ( this->bit(patterntrig[5],0) ){ pmtpattern++;  if ( !this->checkPMT("S11_1A")){ pmtnosignal++; good = false;}}
1648      if ( this->bit(patterntrig[5],1) ){ pmtpattern++;  if ( !this->checkPMT("S11_2A")){ pmtnosignal++; good = false;}}
1649      if ( this->bit(patterntrig[5],2) ){ pmtpattern++;  if ( !this->checkPMT("S11_3A")){ pmtnosignal++; good = false;}}
1650      if ( this->bit(patterntrig[5],3) ){ pmtpattern++;  if ( !this->checkPMT("S11_4A")){ pmtnosignal++; good = false;}}
1651      if ( this->bit(patterntrig[5],4) ){ pmtpattern++;  if ( !this->checkPMT("S11_5A")){ pmtnosignal++; good = false;}}
1652      if ( this->bit(patterntrig[5],5) ){ pmtpattern++;  if ( !this->checkPMT("S11_6A")){ pmtnosignal++; good = false;}}
1653      if ( this->bit(patterntrig[5],6) ){ pmtpattern++;  if ( !this->checkPMT("S11_7A")){ pmtnosignal++; good = false;}}
1654      if ( this->bit(patterntrig[5],7) ){ pmtpattern++;  if ( !this->checkPMT("S11_8A")){ pmtnosignal++; good = false;}}      
1655      if ( this->bit(patterntrig[5],8) ){ pmtpattern++;  if ( !this->checkPMT("S11_1B")){ pmtnosignal++; good = false;}}
1656      if ( this->bit(patterntrig[5],9) ){ pmtpattern++;  if ( !this->checkPMT("S11_2B")){ pmtnosignal++; good = false;}}
1657      if ( this->bit(patterntrig[5],10) ){ pmtpattern++; if ( !this->checkPMT("S11_3B")){ pmtnosignal++; good = false;}}
1658      if ( this->bit(patterntrig[5],11) ){ pmtpattern++; if ( !this->checkPMT("S11_4B")){ pmtnosignal++; good = false;}}
1659      if ( this->bit(patterntrig[5],12) ){ pmtpattern++; if ( !this->checkPMT("S11_5B")){ pmtnosignal++; good = false;}}
1660      if ( this->bit(patterntrig[5],13) ){ pmtpattern++; if ( !this->checkPMT("S11_6B")){ pmtnosignal++; good = false;}}
1661      if ( this->bit(patterntrig[5],14) ){ pmtpattern++; if ( !this->checkPMT("S11_7B")){ pmtnosignal++; good = false;}}
1662      if ( this->bit(patterntrig[5],15) ){ pmtpattern++; if ( !this->checkPMT("S11_8B")){ pmtnosignal++; good = false;}}
1663    
1664      return good;
1665    }
1666    
1667    bool ToFLevel2::checkPMTpmttrig(TrigLevel2 *trg){
1668      //  UInt_t *patterntrig = trg->patterntrig;
1669      int rS11 = 0;
1670      int rS12 = 0;
1671      int rS21 = 0;
1672      int rS22 = 0;
1673      int rS31 = 0;
1674      int rS32 = 0;
1675    
1676      // trigger configuration for the event from saved pmts
1677      TClonesArray* Pmt = this->PMT;
1678      for(int i=0; i<Pmt->GetEntries(); i++) {  
1679        ToFPMT* pmthit = (ToFPMT*)Pmt->At(i);
1680        TString pmtname = this->GetPMTName(pmthit->pmt_id);
1681        if ( pmtname.Contains("S11") ) rS11++;
1682        if ( pmtname.Contains("S12") ) rS12++;
1683        if ( pmtname.Contains("S21") ) rS21++;
1684        if ( pmtname.Contains("S22") ) rS22++;
1685        if ( pmtname.Contains("S31") ) rS31++;
1686        if ( pmtname.Contains("S32") ) rS32++;
1687      }
1688      int rTOF1 = (rS11 + rS12) * (rS21 + rS22) * (rS31 + rS32);
1689      int rTOF2 = (rS11 * rS12) * (rS21 * rS22) * (rS31 * rS32);
1690    
1691      int rTOF3 = (rS21 + rS22) * (rS31 + rS32);
1692      int rTOF4 = (rS21 * rS22) * (rS31 * rS32);
1693    
1694      int rTOF5 = rS12 * (rS21 * rS22);
1695    
1696      int rTOF6 = (rS11 + rS12) * (rS31 + rS32);
1697      int rTOF7 = (rS11 * rS12) * (rS31 * rS32);
1698    
1699    
1700      // trigger configuration of the run
1701      bool TCTOF1 = false;
1702      bool TCTOF2 = false;
1703      bool TCTOF3 = false;
1704      bool TCTOF4 = false;
1705      bool TCTOF5 = false;
1706      bool TCTOF6 = false;
1707      bool TCTOF7 = false;
1708      if ( trg->trigconf & (1<<0) ) TCTOF1 = true;
1709      if ( trg->trigconf & (1<<1) ) TCTOF2 = true;
1710      if ( trg->trigconf & (1<<2) ) TCTOF3 = true;
1711      if ( trg->trigconf & (1<<3) ) TCTOF4 = true;
1712      if ( trg->trigconf & (1<<4) ) TCTOF5 = true;
1713      if ( trg->trigconf & (1<<5) ) TCTOF6 = true;
1714      if ( trg->trigconf & (1<<6) ) TCTOF7 = true;
1715    
1716      // do patterntrig pmts match the trigger configuration?
1717      bool pmtsconf_trigconf_match = true;
1718      if ( rTOF1 == 0 && TCTOF1 ) pmtsconf_trigconf_match = false;
1719      if ( rTOF2 == 0 && TCTOF2 ) pmtsconf_trigconf_match = false;
1720      if ( rTOF3 == 0 && TCTOF3 ) pmtsconf_trigconf_match = false;
1721      if ( rTOF4 == 0 && TCTOF4 ) pmtsconf_trigconf_match = false;
1722      if ( rTOF5 == 0 && TCTOF5 ) pmtsconf_trigconf_match = false;
1723      if ( rTOF6 == 0 && TCTOF6 ) pmtsconf_trigconf_match = false;
1724      if ( rTOF7 == 0 && TCTOF7 ) pmtsconf_trigconf_match = false;
1725    
1726      return pmtsconf_trigconf_match;
1727    }
1728    
1729    void ToFLevel2::printPMT(){
1730      TClonesArray* Pmt = this->PMT;
1731      for(int i=0; i<Pmt->GetEntries(); i++) {  
1732        ToFPMT* pmthit = (ToFPMT*)Pmt->At(i);
1733        TString pmtname = this->GetPMTName(pmthit->pmt_id);
1734        printf(" PMT hit: %s \n",pmtname.Data());
1735      }
1736    }
1737    
1738    
1739    ToFdEdx::ToFdEdx()
1740    {
1741      memset(conn,0,12*sizeof(Bool_t));
1742      memset(ts,0,12*sizeof(UInt_t));
1743      memset(te,0,12*sizeof(UInt_t));
1744      eDEDXpmt = new TArrayF(48);
1745      Define_PMTsat();
1746      Clear();
1747    }
1748    
1749    ToFdEdx::~ToFdEdx(){
1750      Clear();
1751      Delete();
1752    }
1753    
1754    void ToFdEdx::Delete(Option_t *option){
1755      if ( eDEDXpmt ){
1756        eDEDXpmt->Set(0);
1757        if ( eDEDXpmt) delete eDEDXpmt;
1758      }
1759    }
1760    
1761    //------------------------------------------------------------------------
1762    void ToFdEdx::CheckConnectors(UInt_t atime, GL_PARAM *glparam, TSQLServer *dbc)
1763    {
1764      for(int i=0; i<12; i++){
1765        if(atime<=ts[i] || atime>te[i]){
1766          Int_t error=glparam->Query_GL_PARAM(atime,210+i,dbc); // parameters stored in DB in GL_PRAM table
1767          if ( error<0 ) {
1768            conn[i]=false;
1769            ts[i]=0;
1770            te[i]=numeric_limits<UInt_t>::max();
1771          };
1772          if ( !error ){
1773            conn[i]=true;
1774            ts[i]=glparam->FROM_TIME;
1775            te[i]=glparam->TO_TIME;
1776          }
1777          if ( error>0 ){
1778            conn[i]=false;
1779            ts[i]=glparam->TO_TIME;
1780            TSQLResult *pResult;
1781            TSQLRow *row;
1782            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);
1783            pResult=dbc->Query(query.Data());
1784            if(!pResult->GetRowCount()){
1785              te[i]=numeric_limits<UInt_t>::max();
1786            }else{
1787              row=pResult->Next();
1788              te[i]=(UInt_t)atoll(row->GetField(0));
1789            }
1790          }
1791          //
1792          
1793        }
1794      }
1795    
1796    }
1797    //------------------------------------------------------------------------
1798    void ToFdEdx::Clear(Option_t *option)
1799    {
1800      //
1801      // Set arrays and initialize structure
1802      //  eDEDXpmt.Set(48);    eDEDXpmt.Reset(-1);   // Set array size  and reset structure
1803      eDEDXpmt->Set(48);    eDEDXpmt->Reset(-1);   // Set array size  and reset structure
1804      //
1805    };
1806    
1807    //------------------------------------------------------------------------
1808    void ToFdEdx::Print(Option_t *option)
1809    {
1810      //
1811      printf("========================================================================\n");
1812    
1813    };
1814    
1815    //------------------------------------------------------------------------
1816    void ToFdEdx::Init(pamela::tof::TofEvent *tofl0)
1817    {
1818      //
1819      ToFLevel2 tf;
1820      for (Int_t gg=0; gg<4;gg++){
1821        for (Int_t hh=0; hh<12;hh++){
1822          //          tofinput_.tdc[hh][gg]=tofEvent->tdc[gg][hh];          
1823          int mm = tf.GetPMTid(gg,hh);        
1824          adc[mm]= (0xFFF & tofl0->adc[gg][hh]); // EM, exclude warning bits
1825        };      
1826      };
1827      
1828    };
1829    
1830    //------------------------------------------------------------------------
1831    void ToFdEdx::Init(Int_t gg, Int_t hh, Float_t adce)
1832    {
1833      //
1834      ToFLevel2 tf;
1835      //  for (Int_t gg=0; gg<4;gg++){
1836      //    for (Int_t hh=0; hh<12;hh++){
1837      int mm = tf.GetPMTid(gg,hh);    
1838      adc[mm]=adce;
1839      
1840    };
1841    //------------------------------------------------------------------------
1842    void ToFdEdx::Process(UInt_t atime, Float_t betamean, Float_t *xtr_tof, Float_t *ytr_tof, Int_t exitat)
1843    {
1844      bool debug = false;
1845      if ( debug ) printf(" INSIDE TOFDEDX PROCESS \n");
1846      // the parameters should be already initialised by InitPar()
1847      //  printf(" in process \n");
1848      Clear();
1849    
1850     // define angle:  
1851      double dx   = xtr_tof[1] - xtr_tof[5];
1852      double dy   = ytr_tof[0] - ytr_tof[4];
1853      double dr   = sqrt(dx*dx+dy*dy);
1854      double theta=atan(dr/76.81);
1855      //
1856      if ( xtr_tof[1] > 99. ||  xtr_tof[5] > 99. || ytr_tof[0] > 99. ||  ytr_tof[4] > 99. ) theta = 0.;
1857      for (Int_t ii=0; ii<6; ii++){
1858        if ( xtr_tof[ii] > 99. ) xtr_tof[ii] = 0.;
1859        if ( ytr_tof[ii] > 99. ) ytr_tof[ii] = 0.;
1860      };
1861      //
1862      if ( debug ) printf(" theta %f \n",theta);
1863      if ( debug ) printf(" xtr_tof %.1f %.1f %.1f %.1f %.1f %.1f \n",xtr_tof[0],xtr_tof[1],xtr_tof[2],xtr_tof[3],xtr_tof[4],xtr_tof[5]);
1864      if ( debug ) printf(" ytr_tof %.1f %.1f %.1f %.1f %.1f %.1f \n",ytr_tof[0],ytr_tof[1],ytr_tof[2],ytr_tof[3],ytr_tof[4],ytr_tof[5]);
1865      //--------------------- TABLE OF PERIODS WITH HV PROBLEMS ----------------------------
1866      
1867      int Aconn=conn[0];    // PMT 0,20,22,24
1868      int Bconn=conn[1];    // PMT 6,12,26,34
1869      int Cconn=conn[2];    // PMT 4,14,28,32
1870      int Dconn=conn[3];    // PMT 2,8,10,30
1871      int Econn=conn[4];    // PMT 42,43,44,47
1872      int Fconn=conn[5];    // PMT 7,19,23,27
1873      int Gconn=conn[6];    // PMT 3,11,25,33
1874      int Hconn=conn[7];    // PMT 1,9,13,21
1875      int Iconn=conn[8];    // PMT 5,29,31,35
1876      int Lconn=conn[9];    // PMT 37,40,45,46
1877      int Mconn=conn[10];    // PMT 15,16,17,18
1878      int Nconn=conn[11];    // PMT 36,38,39,41
1879      if( false ) cout << Gconn << Iconn << Lconn <<endl; // to avoid compilation warnings
1880        
1881      //  printf(" size %i \n",eDEDXpmt.GetSize());
1882      for( int ii=0; ii<48; ii++ ) {
1883        //
1884        //    eDEDXpmt.SetAt(-1.,ii);
1885        //    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]);
1886        if ( debug ) printf("II %i adc %f \n",ii,adc[ii]);
1887    
1888        if( adc[ii] >= 4095. ){
1889          //      eDEDXpmt[ii] = 0.;
1890          eDEDXpmt->AddAt(0.,ii);
1891          if ( debug ) printf(" %i adc>4095 \n",ii);
1892          continue; // EMILIANO
1893        };
1894    
1895        if( adc[ii] >= (PMTsat[ii]-5.) && adc[ii] < 4095. ){
1896          eDEDXpmt->AddAt(1000.,ii);
1897          if ( debug ) printf(" %i adc> pmtsat && adc<4095 \n",ii);
1898          continue; // EMILIANO
1899        };
1900    
1901        if( adc[ii] <= 0. ) {
1902          eDEDXpmt->AddAt(1500.,ii);
1903          if ( debug ) printf(" %i adc<=0 \n",ii);
1904          continue;
1905        };
1906        //
1907        double adcpC   = f_adcPC( adc[ii] );    // - adc conversion in pC
1908        if ( exitat == 0 ){
1909          eDEDXpmt->AddAt((Float_t)adcpC,ii);
1910          continue;
1911        }
1912        //    printf(" e qua? \n");
1913    
1914        double adccorr = adcpC*fabs(cos(theta));    
1915        if ( debug ) printf(" adccorr %f \n",adccorr);
1916        if(adccorr<=0.){
1917          if ( debug ) printf(" %i adccorr<=0 \n",ii);
1918          //      eDEDXpmt->AddAt((Float_t)adcpC,ii);//?
1919          continue;
1920        }
1921        if ( exitat == 1 ){
1922          eDEDXpmt->AddAt((Float_t)adccorr,ii);
1923          continue;
1924        }
1925        //    printf(" e quo? \n");
1926    
1927        //    int standard=0;
1928        int S115B_ok=0;
1929        int S115B_break=0;
1930    
1931        if(atime<1158720000)S115B_ok=1;
1932        else S115B_break=1;
1933    
1934    
1935        //------------------------------------------------------------------------
1936        //    printf(" e qui? \n");
1937        //---------------------------------------------------- Z reconstruction
1938    
1939        double adcHe, adcnorm, adclin, dEdx;//, Zeta; // EM GCC4.7
1940    
1941        adcHe=-2;
1942        adcnorm=-2;
1943        adclin=-2;
1944        dEdx=-2;
1945        //    Zeta=-2;//EM GCC4.7
1946        Double_t correction = 1.;
1947    
1948        if(Aconn==1 && (ii==0 || ii==20 || ii==22 || ii==24)){
1949          correction = 1.675;
1950        }
1951        else if(Bconn==1 && (ii==6 || ii==12 || ii==26 || ii==34)){
1952          correction = 2.482;
1953        }
1954        else if(Cconn==1 && (ii==4 || ii==14 || ii==28 || ii==32)){
1955          correction = 1.464;
1956        }
1957        else if(Dconn==1 && (ii==2 || ii==8 || ii==10 || ii==30)){
1958          correction = 1.995;
1959        }
1960        else if(Econn==1 && (ii==42 || ii==43 || ii==44 || ii==47)){
1961          correction = 1.273;
1962        }
1963        else if(Fconn==1 && (ii==7 || ii==19 || ii==23 || ii==27)){
1964          correction = 1.565;
1965        }
1966        else if(Mconn==1 && (ii==15 || ii==16 || ii==17 || ii==18)){
1967          correction = 1.565;
1968        }
1969        else if(Nconn==1 && (ii==36 || ii==38 || ii==39 || ii==41)){
1970          correction = 1.018;
1971        }
1972        else if(Hconn==1 && (ii==1 || ii==13 || ii==21 || (ii==9&&S115B_ok==1))){
1973          correction = 1.84;
1974        }
1975        else if(S115B_break==1 && ii==9 && Hconn==1){
1976          correction = 1.64;
1977        }
1978        else correction = 1.;
1979        
1980        if( ii==9 && S115B_break==1 ){
1981          adcHe   = f_att5B( ytr_tof[0] )/correction;
1982        } else {
1983          adcHe   = Get_adc_he(ii, xtr_tof, ytr_tof)/correction;
1984        };
1985        if(adcHe<=0){
1986          if ( debug ) printf(" %i adcHe<=0 \n",ii);
1987          //      eDEDXpmt->AddAt((Float_t)adccorr,ii); //?
1988          continue;
1989        }
1990        if ( exitat == 2 ){
1991          if(ii==9 && S115B_break==1)  eDEDXpmt->AddAt(36.*(Float_t)adccorr/adcHe,ii);
1992          else  adclin  = 4.*(Float_t)adccorr/adcHe;
1993          continue;
1994        }
1995    
1996        if(ii==9 && S115B_break==1)  adcnorm = f_pos5B(adccorr);
1997        else adcnorm = f_pos( (parPos[ii]), adccorr);
1998        if(adcnorm<=0){
1999          if ( debug ) printf(" %i adcnorm<=0 \n",ii);
2000          //      eDEDXpmt->AddAt((Float_t)adccorr,ii);//?
2001          continue;
2002        }
2003        if ( debug ) printf(" adcnorm %f \n",adcnorm);
2004    
2005        if(ii==9 && S115B_break==1)  adclin  = 36.*adcnorm/adcHe;
2006        else  adclin  = 4.*adcnorm/adcHe;
2007        if ( debug ) printf(" adclin %f \n",adclin);
2008        if(adclin<=0){
2009          if ( debug ) printf(" %i adclin<=0 \n",ii);
2010          //      eDEDXpmt->AddAt((Float_t)adccorr,ii);//?
2011          continue;
2012        }
2013        if ( exitat == 3 ){
2014          if(ii==9 && S115B_break==1)  eDEDXpmt->AddAt((Float_t)adclin,ii);
2015          else  eDEDXpmt->AddAt((Float_t)adclin,ii);
2016          continue;
2017        }
2018        //
2019        if ( betamean > 99. ){
2020          //      eDEDXpmt.AddAt((Float_t)adclin,ii);
2021          eDEDXpmt->AddAt((Float_t)adclin,ii);
2022          //      printf(" AAPMT IS %i dedx is %f vector is %f \n",ii,adclin,eDEDXpmt[ii]);
2023          if ( debug ) printf(" %i betamean > 99 \n",ii);
2024          continue;
2025        };
2026        //
2027        double dEdxHe=-2;
2028        if(ii==9 && S115B_break==1){
2029          if( betamean <1. ) dEdxHe = f_BB5B( betamean );
2030          else                       dEdxHe = 33;
2031        } else {
2032          if( betamean <1. ) dEdxHe = f_BB( (parBBneg[ii]), betamean );
2033          else                       dEdxHe = parBBpos[ii];
2034        }
2035        
2036        if ( debug ) printf(" dEdxHe %f \n",dEdxHe);
2037        
2038        if(dEdxHe<=0){
2039          eDEDXpmt->AddAt((Float_t)adclin,ii);
2040          if ( debug ) printf(" %i dEdxHe<=0 \n",ii);
2041          continue;
2042        };
2043    
2044        if(ii==9 && S115B_break==1)  dEdx = f_desatBB5B( adclin );
2045        else  dEdx = f_desatBB((parDesatBB[ii]), adclin );
2046    
2047        if(dEdx<=0){
2048          eDEDXpmt->AddAt((Float_t)adclin,ii);
2049          if ( debug ) printf(" %i dEdx<=0 \n",ii);
2050          continue;
2051        };
2052    
2053        if ( debug ) printf(" dEdx %f \n",dEdx);
2054        eDEDXpmt->AddAt((Float_t)dEdx,ii);
2055        //    eDEDXpmt.AddAt((Float_t)dEdx,ii);
2056    
2057        //    printf(" PMT IS %i dedx is %f vector is %f \n",ii,dEdx,eDEDXpmt[ii]);
2058    
2059      }  //end loop on 48 PMT
2060    
2061    };
2062    
2063    
2064    //------------------------------------------------------------------------
2065    void ToFdEdx::Define_PMTsat()
2066    {
2067      Float_t  sat[48] = {
2068        3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
2069        3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
2070        3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
2071        3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
2072        3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
2073        3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };
2074      PMTsat.Set(48,sat);
2075    }
2076    
2077    //------------------------------------------------------------------------
2078    void ToFdEdx::ReadParBBpos( const char *fname )
2079    {
2080      //  printf("read %s\n",fname);
2081      parBBpos.Set(48);
2082      FILE *fattin = fopen( fname , "r" );
2083      for (int i=0; i<48; i++) {
2084        int   tid=0;
2085        float  tp;
2086        if(fscanf(fattin,"%d %f",
2087                  &tid, &tp )!=2) break;
2088        parBBpos[i]=tp;
2089      }
2090      fclose(fattin);
2091    }
2092    
2093    //------------------------------------------------------------------------
2094    void ToFdEdx::ReadParDesatBB( const char *fname )
2095    {
2096      //  printf("read %s\n",fname);
2097      FILE *fattin = fopen( fname , "r" );
2098      for (int i=0; i<48; i++) {
2099        int   tid=0;
2100        float  tp[3];
2101        if(fscanf(fattin,"%d %f %f %f",
2102                  &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
2103        parDesatBB[i].Set(3,tp);
2104      }
2105      fclose(fattin);
2106    }
2107    
2108    
2109    //------------------------------------------------------------------------
2110    void ToFdEdx::ReadParBBneg( const char *fname )
2111    
2112    {
2113      //  printf("read %s\n",fname);
2114      FILE *fattin = fopen( fname , "r" );
2115      for (int i=0; i<48; i++) {
2116        int   tid=0;
2117        float  tp[3];
2118        if(fscanf(fattin,"%d %f %f %f",
2119                  &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
2120        parBBneg[i].Set(3,tp);
2121      }
2122      fclose(fattin);
2123    }
2124    
2125    //------------------------------------------------------------------------
2126    void ToFdEdx::ReadParPos( const char *fname )
2127    {
2128      //  printf("read %s\n",fname);
2129      FILE *fattin = fopen( fname , "r" );
2130      for (int i=0; i<48; i++) {
2131        int   tid=0;
2132        float  tp[4];
2133        if(fscanf(fattin,"%d %f %f %f %f",
2134                  &tid, &tp[0], &tp[1], &tp[2], &tp[3])!=5) break;
2135        parPos[i].Set(4,tp);
2136      }
2137      fclose(fattin);
2138    }
2139    
2140    //------------------------------------------------------------------------
2141    void ToFdEdx::ReadParAtt( const char *fname )
2142    {
2143      //  printf("read %s\n",fname);
2144      FILE *fattin = fopen( fname , "r" );
2145      for (int i=0; i<48; i++) {
2146        int   tid=0;
2147        float  tp[6];
2148        if(fscanf(fattin,"%d %f %f %f %f %f %f",
2149                  &tid, &tp[0], &tp[1], &tp[2], &tp[3], &tp[4], &tp[5] )!=7) break;
2150        parAtt[i].Set(6,tp);
2151      }
2152      fclose(fattin);
2153    }
2154    
2155    
2156    
2157    
2158    
2159    
2160    double ToFdEdx::f_att( TArrayF &p, float x )
2161    {
2162      return
2163        p[0] +
2164        p[1]*x +
2165        p[2]*x*x +
2166        p[3]*x*x*x +
2167        p[4]*x*x*x*x +
2168        p[5]*x*x*x*x*x;
2169    }
2170    //------------------------------------------------------------------------
2171    double ToFdEdx::f_att5B( float x )
2172    {
2173      return
2174        101.9409 +
2175        6.643781*x +
2176        0.2765518*x*x +
2177        0.004617647*x*x*x +
2178        0.0006195132*x*x*x*x +
2179        0.00002813734*x*x*x*x*x;
2180    }
2181    
2182    
2183    double ToFdEdx::f_pos( TArrayF &p, float x )
2184    {
2185      return
2186        p[0] +
2187        p[1]*x +
2188        p[2]*x*x +
2189        p[3]*x*x*x;
2190    }
2191    
2192    double ToFdEdx::f_pos5B( float x )
2193    {
2194      return
2195        15.45132 +
2196        0.8369721*x +
2197        0.0005*x*x;
2198    }
2199    
2200    
2201    
2202    double ToFdEdx::f_adcPC( float x )
2203    {
2204      return 28.12+0.6312*x-5.647e-05*x*x+3.064e-08*x*x*x;
2205    }
2206    
2207    
2208    float ToFdEdx::Get_adc_he( int id, float pl_x[6], float pl_y[6])
2209    {
2210    
2211      //
2212      // input: id - pmt [0:47}
2213      //             pl_x - coord x of the tof plane
2214      //             pl_y - coord y
2215    
2216      adc_he = 0;
2217      if( eGeom.GetXY(id)==1 )  adc_he = f_att( (parAtt[id]), pl_x[eGeom.GetPlane(id)] );
2218      if( eGeom.GetXY(id)==2 )  adc_he = f_att( (parAtt[id]), pl_y[eGeom.GetPlane(id)] );
2219      return adc_he;
2220    }
2221    
2222    //------------------------------------------------------------------------
2223    double ToFdEdx::f_BB( TArrayF &p, float x )
2224    {
2225      return  p[0]/(x*x)*(log(x*x/(1-x*x)) - p[1]*x*x - p[2]);
2226    }
2227    
2228    //------------------------------------------------------------------------
2229    double ToFdEdx::f_BB5B( float x )
2230    {
2231      return  0.165797/(x*x)*(log(x*x/(1-x*x)) + 140.481*x*x + 52.9258);
2232    }
2233    //------------------------------------------------------------------------
2234    double ToFdEdx::f_desatBB( TArrayF &p, float x )
2235    {
2236      return
2237        p[0] +
2238        p[1]*x +
2239        p[2]*x*x;
2240    }
2241    
2242    //------------------------------------------------------------------------
2243    double ToFdEdx::f_desatBB5B( float x )
2244    {
2245      return
2246        -2.4 +
2247        0.75*x +
2248        0.009*x*x;
2249    }
2250    

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.47

  ViewVC Help
Powered by ViewVC 1.1.23