/** * \file ToFLevel2.cpp * \author Gianfranca DeRosa, Wolfgang Menn * * WM dec 2008: Description of "GetdEdx" changed * WM dec 2008: "GetdEdxPaddle" modified: Now includes saturation limit * PMTs higher than the saturation limit are not used for dEdx * WM apr 2009: bug found by Nicola in method "GetPaddlePlane" */ #include using namespace std; ClassImp(ToFPMT); ClassImp(ToFdEdx); ClassImp(ToFGeom); ClassImp(ToFTrkVar); ClassImp(ToFLevel2); ToFPMT::ToFPMT(){ pmt_id = 0; adc = 0.; tdc_tw = 0.; tdc = 0.; } ToFPMT::ToFPMT(const ToFPMT &t){ pmt_id = t.pmt_id; adc = t.adc; tdc_tw = t.tdc_tw; tdc = t.tdc; } void ToFPMT::Clear(Option_t *t){ pmt_id = 0; adc = 0.; tdc_tw = 0.; tdc = 0.; } ToFTrkVar::ToFTrkVar() { trkseqno = 0; npmttdc = 0; npmtadc = 0; pmttdc = TArrayI(48); pmtadc = TArrayI(48); tdcflag = TArrayI(48); // gf: 30 Nov 2006 adcflag = TArrayI(48); // gf: 30 Nov 2006 dedx = TArrayF(48); // // memset(beta, 0, 13*sizeof(Float_t)); memset(xtofpos, 0, 3*sizeof(Float_t)); memset(ytofpos, 0, 3*sizeof(Float_t)); memset(xtr_tof, 0, 6*sizeof(Float_t)); memset(ytr_tof, 0, 6*sizeof(Float_t)); // }; void ToFTrkVar::Clear(Option_t *t) { trkseqno = 0; npmttdc = 0; npmtadc = 0; pmttdc.Reset(); pmtadc.Reset(); tdcflag.Reset(); // gf: 30 Nov 2006 adcflag.Reset(); // gf: 30 Nov 2006 dedx.Reset(); // memset(beta, 0, 13*sizeof(Float_t)); memset(xtofpos, 0, 3*sizeof(Float_t)); memset(ytofpos, 0, 3*sizeof(Float_t)); memset(xtr_tof, 0, 6*sizeof(Float_t)); memset(ytr_tof, 0, 6*sizeof(Float_t)); // }; ToFTrkVar::ToFTrkVar(const ToFTrkVar &t){ trkseqno = t.trkseqno; // npmttdc = t.npmttdc; npmtadc = t.npmtadc; (t.pmttdc).Copy(pmttdc); (t.pmtadc).Copy(pmtadc); (t.tdcflag).Copy(tdcflag); // gf: 30 Nov 2006 (t.adcflag).Copy(adcflag); // gf: 30 Nov 2006 (t.dedx).Copy(dedx); // memcpy(beta,t.beta,sizeof(beta)); memcpy(xtofpos,t.xtofpos,sizeof(xtofpos)); memcpy(ytofpos,t.ytofpos,sizeof(ytofpos)); memcpy(xtr_tof,t.xtr_tof,sizeof(xtr_tof)); memcpy(ytr_tof,t.ytr_tof,sizeof(ytr_tof)); // }; ToFLevel2::ToFLevel2() { // // PMT = new TClonesArray("ToFPMT",12); //ELENA // ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA PMT = 0; //ELENA ToFTrk = 0; //ELENA // this->Clear(); // }; void ToFLevel2::Set(){//ELENA if(!PMT)PMT = new TClonesArray("ToFPMT",12); //ELENA if(!ToFTrk)ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA }//ELENA void ToFLevel2::Clear(Option_t *t){ // if(ToFTrk)ToFTrk->Delete(); //ELENA if(PMT)PMT->Delete(); //ELENA memset(tof_j_flag, 0, 6*sizeof(Int_t)); unpackError = 0; // }; void ToFLevel2::Delete(Option_t *t){ //ELENA // if(ToFTrk){ ToFTrk->Delete(); //ELENA delete ToFTrk; //ELENA } if(PMT){ PMT->Delete(); //ELENA delete PMT; //ELENA } //ELENA // }; //ELENA ToFTrkVar *ToFLevel2::GetToFTrkVar(Int_t itrk){ // if(itrk >= ntrk()){ printf(" ToFLevel2 ERROR: track related variables set %i does not exists! \n",itrk); printf(" stored track related variables = %i \n",ntrk()); return(NULL); } // if(!ToFTrk)return 0; //ELENA TClonesArray &t = *(ToFTrk); ToFTrkVar *toftrack = (ToFTrkVar*)t[itrk]; return toftrack; } ToFPMT *ToFLevel2::GetToFPMT(Int_t ihit){ // if(ihit >= npmt()){ printf(" ToFLevel2 ERROR: pmt variables set %i does not exists! \n",ihit); printf(" stored pmt variables = %i \n",npmt()); return(NULL); } // if(!PMT)return 0; //ELENA TClonesArray &t = *(PMT); ToFPMT *tofpmt = (ToFPMT*)t[ihit]; return tofpmt; } //-------------------------------------- // // //-------------------------------------- /** * Method to get the plane ID (11 12 21 22 31 32) from the plane index (0 1 2 3 4 5) * @param Plane index (0,1,2,3,4,5). */ Int_t ToFLevel2::GetToFPlaneID(Int_t ip){ if(ip>=0 && ip<6)return 10*((int)(ip/2+1.1))+(ip%2)+1; else return -1; }; /** * Method to get the plane index (0 1 2 3 4 5) from the plane ID (11 12 21 22 31 32) * @param plane Plane ID (11, 12, 21, 22, 31, 32) */ Int_t ToFLevel2::GetToFPlaneIndex(Int_t plane_id){ if( plane_id == 11 || plane_id == 12 || plane_id == 21 || plane_id == 22 || plane_id == 31 || plane_id == 32 || false)return (Int_t)(plane_id/10)*2-1- plane_id%2; else return -1; }; /** * Method to know if a given ToF paddle was hit, that is there is a TDC signal * from both PMTs. The method uses the "tof_j_flag" variable. * @param plane Plane ID (11, 12, 21, 22, 31, 32) or Plane index (0,1,2,3,4,5). * @param paddle_id Paddle ID. * @return 1 if the paddle was hit. */ Bool_t ToFLevel2::HitPaddle(Int_t plane, Int_t paddle_id){ //<<< NEW Int_t ip = -1; if (plane>=6 ) ip = GetToFPlaneIndex(plane); else if(plane>=0 && plane < 6) ip = plane; Int_t flag=0; if(ip != -1)flag = tof_j_flag[ip] & (int)pow(2.,(double)paddle_id); if( (ip == 0 && paddle_id < 8 && flag) || (ip == 1 && paddle_id < 6 && flag) || (ip == 2 && paddle_id < 2 && flag) || (ip == 3 && paddle_id < 2 && flag) || (ip == 4 && paddle_id < 3 && flag) || (ip == 5 && paddle_id < 3 && flag) || false) return true; else return false; }; /** * Method to get the number of hit paddles on a ToF plane. * @param plane Plane ID (11, 12, 21, 22, 31, 32) or Plane index (0,1,2,3,4,5). */ Int_t ToFLevel2::GetNHitPaddles(Int_t plane){ Int_t npad=0; for(Int_t i=0; i<8; i++)npad = npad + (int)HitPaddle(plane,i); return npad; }; //wm Nov 08 //gf Apr 07 /** * Method to get the mean dEdx from a ToF layer - ATTENTION: * It will sum up the dEdx of all the paddles, but since by definition * only the paddle hitted by the track gets a dEdx value and the other * paddles are set to zero, the output is just the dEdx of the hitted * paddle in each layer! * The "adcfl" option is not very useful (an artificial dEdx is per * definition= 1 mip and not a real measurement), anyway left in the code * @param notrack Track Number * @param plane Plane index (0,1,2,3,4,5) * @param adcflag in the plane (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; ) */ Float_t ToFLevel2::GetdEdx(Int_t notrack, Int_t plane, Int_t adcfl){ Float_t dedx = 0.; Float_t PadEdx =0.; Int_t SatWarning; Int_t pad=-1; // ToFTrkVar *trk = GetToFTrkVar(notrack); if(!trk) return 0; //ELENA // for (Int_t ii=0; iinpmtadc; i++){ // pmt_id = (trk->pmtadc).At(i); // GetPMTIndex(pmt_id,hh,kk); adc[kk][hh] = (trk->dedx).At(i); // }; // for (Int_t i=0; ipmt_id,hh,kk); // tdc[kk][hh] = pmt->tdc_tw; // }; // return; }; /** * Method to get the plane index (0 - 5) for the PMT_ID as input * @param pmt_id PMT_ID (0 - 47) */ Int_t ToFLevel2::GetPlaneIndex(Int_t pmt_id){ TString pmtname = GetPMTName(pmt_id); pmtname.Resize(3); if ( !strcmp(pmtname,"S11") ) return(0); if ( !strcmp(pmtname,"S12") ) return(1); if ( !strcmp(pmtname,"S21") ) return(2); if ( !strcmp(pmtname,"S22") ) return(3); if ( !strcmp(pmtname,"S31") ) return(4); if ( !strcmp(pmtname,"S32") ) return(5); return(-1); }; /** * Method to get the PMT_ID if the index (4,12) is given. We have 4 channels on * each of the 12 half-boards, this method decodes which PMT is cables to which * channel. * @param hh Channel * @param kk HalfBoard */ Int_t ToFLevel2::GetPMTid(Int_t hh, Int_t kk){ // short tof[4][24] = { {4, 4, 4, 4, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 2, 3, 3, 3, 3, 4}, {1, 3, 5, 7, 10, 12, 2, 4, 2, 4, 6, 8, 10, 12, 1, 5, 3, 9, 7, 9, 11, 1, 5, 9}, {2, 2, 2, 2, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 2, 1, 2, 1, 2, 2, 2, 3, 3, 4}, {6, 8, 12, 10, 8, 6, 4, 2, 12, 10, 8, 6, 4, 2, 9, 7, 11, 11, 5, 3, 1, 3, 7, 11} }; // Int_t ind = 0; Int_t k = 0; while (k < 24){ Int_t j = 0; while (j < 2){ Int_t ch = tof[2*j][k] - 1; Int_t hb = tof[2*j + 1][k] - 1; /* tofEvent->tdc[ch][hb] */ if( ch == hh && hb == kk ){ ind = 2*k + j; break; }; j++; }; k++; }; return ind; }; /** * Method to get the PMT index if the PMT ID is given. This method is the * "reverse" of method "GetPMTid" * @param ind PMT_ID (0 - 47) * @param hb HalfBoard * @param ch Channel */ void ToFLevel2::GetPMTIndex(Int_t ind, Int_t &hb, Int_t &ch){ // short tof[4][24] = { {4, 4, 4, 4, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 2, 3, 3, 3, 3, 4}, {1, 3, 5, 7, 10, 12, 2, 4, 2, 4, 6, 8, 10, 12, 1, 5, 3, 9, 7, 9, 11, 1, 5, 9}, {2, 2, 2, 2, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 2, 1, 2, 1, 2, 2, 2, 3, 3, 4}, {6, 8, 12, 10, 8, 6, 4, 2, 12, 10, 8, 6, 4, 2, 9, 7, 11, 11, 5, 3, 1, 3, 7, 11} }; // Int_t k = 0; while (k < 24){ Int_t j = 0; while (j < 2){ /* tofEvent->tdc[ch][hb] */ if( ind == 2*k + j ){ ch = tof[2*j][k] - 1; hb = tof[2*j + 1][k] - 1; return; }; j++; }; k++; }; return; }; // wm Nov 08 revision - saturation values included /// gf Apr 07 /** * Method to get the dEdx from a given ToF paddle. * If two PMTs are good, the mean dEdx of both PMTs is taken, otherwise * just the dEdx of the "good" PMT. If both PMTs are above saturation => dEdx=1000 * @param notrack Track Number * @param Paddle index (0,1,...,23). * @param adcflag in the paddle (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; ) * @param PadEdx dEdx from a given ToF paddle * @param SatWarning 1 if the PMT ios near saturation region (adcraw ~3000) */ void ToFLevel2::GetdEdxPaddle(Int_t notrack, Int_t paddleid, Int_t adcfl, Float_t &PadEdx, Int_t &SatWarning){ /* Float_t PMTsat[48] = { 3162.14, 3165.48, 3153.85, 3085.73, 3089.65, 3107.64, 3097.52, 3078.37, 3130.05, 3087.07, 3112.22, 3102.92, 3080.58, 3092.55, 3087.94, 3125.03, 3094.09, 3143.16, 3125.51, 3181.27, 3092.09, 3124.98, 3069.3, 3095.53, 3097.11, 3133.53, 3114.73, 3113.01, 3091.19, 3097.99, 3033.84, 3134.98, 3081.37, 3111.04, 3066.77, 3108.17, 3133, 3111.06, 3052.52, 3140.66, 3106.33, 3094.85, 3150.85, 3118.8, 3096.24, 3118.47,3111.36, 3117.11 } ; */ // new values from Napoli dec 2008 Float_t PMTsat[48] = { 3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27, 3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14, 3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62, 3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26, 3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18, 3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 }; for (Int_t i=0; i<48;i++) PMTsat[i] = PMTsat[i] - 5.; // safety margin PadEdx = 0.; // SatWarning = 1000; SatWarning = 0; // 0=good, increase for each bad PMT Float_t dEdx[48] = {0}; Int_t pmt_id = -1; Float_t adcraw[48]; // ToFTrkVar *trk = GetToFTrkVar(notrack); if(!trk) return; //ELENA // Int_t pmtleft=-1; Int_t pmtright=-1; GetPaddlePMT(paddleid, pmtleft, pmtright); adcraw[pmtleft] = 4095; adcraw[pmtright] = 4095; for (Int_t jj=0; jjpmt_id; if(pmt_id==pmtleft){ adcraw[pmtleft] = pmt->adc; } if(pmt_id==pmtright){ adcraw[pmtright] = pmt->adc; } } for (Int_t i=0; inpmtadc; i++){ if((trk->adcflag).At(i)==0 || adcfl==100){ if((trk->pmtadc).At(i) == pmtleft)dEdx[pmtleft] = (trk->dedx).At(i); if((trk->pmtadc).At(i) == pmtright)dEdx[pmtright] = (trk->dedx).At(i); }else{ if((trk->pmtadc).At(i) == pmtleft)dEdx[pmtleft] = 0.; if((trk->pmtadc).At(i) == pmtright)dEdx[pmtright] = 0.; } } // if( adcraw[pmtleft] >3000 || adcraw[pmtright] >3000)SatWarning=1; //old version // Increase SatWarning Counter for each PMT>Sat if( adcraw[pmtleft] > PMTsat[pmtleft])SatWarning++; if( adcraw[pmtright] > PMTsat[pmtright])SatWarning++; // if ADC > sat set dEdx=1000 if( adcraw[pmtleft] > PMTsat[pmtleft]) dEdx[pmtleft] = 1000.; if( adcraw[pmtright] > PMTsat[pmtright]) dEdx[pmtright] = 1000. ; // if two PMT are good, take mean dEdx, otherwise only the good dEdx if(dEdx[pmtleft]<1000 && dEdx[pmtright]<1000) PadEdx = (dEdx[pmtleft]+dEdx[pmtright])*0.5; if(dEdx[pmtleft]==1000 && dEdx[pmtright]<1000) PadEdx = dEdx[pmtright]; if(dEdx[pmtleft]<1000 && dEdx[pmtright]==1000) PadEdx = dEdx[pmtleft]; }; // // gf Apr 07 /** * Method to get the PMT name (like "S11_1A") if the PMT_ID is given. * Indexes of corresponding plane, paddle and pmt are also given as output. * @param ind PMT_ID (0 - 47) * @param iplane plane index (0 - 5) * @param ipaddle paddle index (relative to the plane) * @param ipmt pmt index (0(A), 1(B)) */ TString ToFLevel2::GetPMTName(Int_t ind, Int_t &iplane, Int_t &ipaddle,Int_t &ipmt){ TString pmtname = " "; TString photoS[48] = { "S11_1A", "S11_1B", "S11_2A", "S11_2B", "S11_3A", "S11_3B", "S11_4A", "S11_4B", "S11_5A", "S11_5B", "S11_6A", "S11_6B", "S11_7A", "S11_7B", "S11_8A", "S11_8B", "S12_1A", "S12_1B", "S12_2A", "S12_2B", "S12_3A", "S12_3B", "S12_4A", "S12_4B", "S12_5A", "S12_5B", "S12_6A", "S12_6B", "S21_1A", "S21_1B", "S21_2A", "S21_2B", "S22_1A", "S22_1B", "S22_2A", "S22_2B", "S31_1A", "S31_1B", "S31_2A", "S31_2B", "S31_3A", "S31_3B", "S32_1A", "S32_1B", "S32_2A", "S32_2B", "S32_3A", "S32_3B" }; pmtname = photoS[ind].Data(); TString ss = pmtname(1,2); iplane = (int)(atoi(ss.Data())/10)*2-3+atoi(ss.Data())%10; ss = pmtname(4); ipaddle = atoi(ss.Data())-1 ; if( pmtname.Contains("A") )ipmt=0; if( pmtname.Contains("B") )ipmt=1; return pmtname; }; /** * Method to get the PMT name (like "S11_1A") if the PMT_ID is given * @param ind PMT_ID (0 - 47) */ TString ToFLevel2::GetPMTName(Int_t ind){ Int_t iplane = -1; Int_t ipaddle = -1; Int_t ipmt = -1; return GetPMTName(ind,iplane,ipaddle,ipmt); }; // wm jun 08 Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane){ return GetPaddleIdOfTrack(xtr ,ytr ,plane, 0.4); } // gf Apr 07 Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane, Float_t margin){ Double_t xt,yt,xl,xh,yl,yh; Float_t tof11_x[8] = {-17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85}; Float_t tof12_y[6] = { -13.75,-8.25,-2.75,2.75,8.25,13.75}; Float_t tof21_y[2] = { 3.75,-3.75}; Float_t tof22_x[2] = { -4.5,4.5}; Float_t tof31_x[3] = { -6.0,0.,6.0}; Float_t tof32_y[3] = { -5.0,0.0,5.0}; // S11 8 paddles 33.0 x 5.1 cm // S12 6 paddles 40.8 x 5.5 cm // S21 2 paddles 18.0 x 7.5 cm // S22 2 paddles 15.0 x 9.0 cm // S31 3 paddles 15.0 x 6.0 cm // S32 3 paddles 18.0 x 5.0 cm Int_t paddleidoftrack=-1; // //--- S11 ------ if(plane==0){ xt = xtr; yt = ytr; paddleidoftrack=-1; yl = -33.0/2. ; yh = 33.0/2. ; if ((yt>yl)&&(ytxl)&&(xtxl)&&(xtyl)&&(ytxl)&&(xtyl)&&(ytyl)&&(ytxl)&&(xtyl)&&(ytxl)&&(xtxl)&&(xtyl)&&(yt "x"-sigma. A chi2 value is * calculated, furthermore a "quality" value by adding the weights which * are finally used. If all measurements are taken, "quality" will be = 22.47. * A chi2 cut around 3-4 and a quality-cut > 20 is needed for clean beta * measurements like antiprotons etc. * The Level2 output is derived in the fortran routines using: 10.,10.,20. * @param notrack Track Number * @param cut on residual: difference between single measurement and mean * @param cut on "quality" * @param cut on chi2 */ Float_t ToFLevel2::CalcBeta(Int_t notrack, Float_t resmax, Float_t qualitycut, Float_t chi2cut){ // cout<<" in CalcBeta "<beta[i]; } //======================================================================== //--- Find out ToF layers with artificial TDC values & fill vector --- //======================================================================== Float_t w_il[6]; for (Int_t jj=0; jj<6;jj++) { w_il[jj] = 1000.; } for (Int_t i=0; inpmttdc; i++){ // pmt_id = (trk->pmttdc).At(i); pmt_plane = GetPlaneIndex(pmt_id); tdcfl = (trk->tdcflag).At(i); if (w_il[pmt_plane] != 1.) w_il[pmt_plane] = tdcfl; //tdcflag }; //======================================================================== //--- Set weights for the 12 measurements using information for top and bottom: //--- if no measurements: weight = set to very high value=> not used //--- top or bottom artificial: weight*sqrt(2) //--- top and bottom artificial: weight*sqrt(2)*sqrt(2) //======================================================================== Int_t itop[12] = {0,0,1,1,2,2,3,3,0,0,1,1}; Int_t ibot[12] = {4,5,4,5,4,5,4,5,2,3,2,3}; xhelp= 1E09; for (Int_t jj=0; jj<12;jj++) { if (jj<4) xhelp = 0.11; // S1-S3 if ((jj>3)&&(jj<8)) xhelp = 0.18; // S2-S3 if (jj>7) xhelp = 0.28; // S1-S2 if ((w_il[itop[jj]] == 1000.) && (w_il[ibot[jj]] == 1000.)) xhelp = 1E09; if ((w_il[itop[jj]] == 1) || (w_il[ibot[jj]] == 1.)) xhelp = xhelp*1.414 ; if ((w_il[itop[jj]] == 1) && (w_il[ibot[jj]] == 1.)) xhelp = xhelp*2. ; w_i[jj] = 1./xhelp; } //======================================================================== //--- Calculate mean beta for the first time ----------------------------- //--- We are using "1/beta" since its error is gaussian ------------------ //======================================================================== Int_t icount=0; sw=0.; sxw=0.; beta_mean=100.; for (Int_t jj=0; jj<12;jj++){ if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.)) { icount= icount+1; sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ; sw =sw + w_i[jj]*w_i[jj] ; } } if (icount>0) beta_mean=1./(sxw/sw); beta_mean_inv = 1./beta_mean; //======================================================================== //--- Calculate beta for the second time, use residuals of the single //--- measurements to get a chi2 value //======================================================================== icount=0; sw=0.; sxw=0.; betachi = 100.; chi2 = 0.; quality=0.; for (Int_t jj=0; jj<12;jj++){ if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.)&&(w_i[jj]>0.01)) { res = beta_mean_inv - (1./b[jj]) ; if (fabs(res*w_i[jj])0) chi2 = chi2/(icount) ; if (icount>0) betachi=1./(sxw/sw); bxx = 100.; if ((chi2 < chi2cut)&&(quality>qualitycut)) bxx = betachi; // return(bxx); }; //////////////////////////////////////////////////// //////////////////////////////////////////////////// /** * Fills a struct cToFLevel2 with values from a ToFLevel2 object (to put data into a F77 common). */ void ToFLevel2::GetLevel2Struct(cToFLevel2 *l2) const{ for(Int_t i=0;i<6;i++) l2->tof_j_flag[i]=tof_j_flag[i]; if(ToFTrk){ //ELENA l2->ntoftrk = ToFTrk->GetEntries(); for(Int_t j=0;jntoftrk;j++){ l2->toftrkseqno[j]= ((ToFTrkVar*)ToFTrk->At(j))->trkseqno; l2->npmttdc[j]= ((ToFTrkVar*)ToFTrk->At(j))->npmttdc; for(Int_t i=0;inpmttdc[j];i++){ l2->pmttdc[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->pmttdc.At(i); l2->tdcflag[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->tdcflag.At(i); // gf: 30 Nov 2006 } for(Int_t i=0;i<13;i++) l2->beta[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->beta[i]; l2->npmtadc[j]= ((ToFTrkVar*)ToFTrk->At(j))->npmtadc; for(Int_t i=0;inpmtadc[j];i++){ l2->pmtadc[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->pmtadc.At(i); l2->adcflag[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->adcflag.At(i); // gf: 30 Nov 2006 l2->dedx[i][j] = ((ToFTrkVar*)ToFTrk->At(j))->dedx.At(i); } for(Int_t i=0;i<3;i++){ l2->xtofpos[i][j]=((ToFTrkVar*)ToFTrk->At(j))->xtofpos[i]; l2->ytofpos[i][j]=((ToFTrkVar*)ToFTrk->At(j))->ytofpos[i]; } for(Int_t i=0;i<6;i++){ l2->xtr_tof[i][j]=((ToFTrkVar*)ToFTrk->At(j))->xtr_tof[i]; l2->ytr_tof[i][j]=((ToFTrkVar*)ToFTrk->At(j))->ytr_tof[i]; } } } //ELENA if(PMT){ //ELENA l2->npmt = PMT->GetEntries(); for(Int_t j=0;jnpmt;j++){ l2->pmt_id[j] = ((ToFPMT*)PMT->At(j))->pmt_id; l2->adc[j] =((ToFPMT*)PMT->At(j))->adc; l2->tdc_tw[j] =((ToFPMT*)PMT->At(j))->tdc_tw; } } //ELENA } // // Reprocessing tool // Emiliano 08/04/07 // Int_t ToFLevel2::Process(TrkLevel2 *trk, TrigLevel2 *trg, GL_RUN *run, OrbitalInfo *orb, Bool_t force){ // // Copiare qui qualcosa di simile a calonuclei per evitare di riprocessare sempre tutto // // // structures to communicate with F77 // extern struct ToFInput tofinput_; extern struct ToFOutput tofoutput_; // // DB connection // TString host; TString user; TString psw; const char *pamdbhost=gSystem->Getenv("PAM_DBHOST"); const char *pamdbuser=gSystem->Getenv("PAM_DBUSER"); const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW"); if ( !pamdbhost ) pamdbhost = ""; if ( !pamdbuser ) pamdbuser = ""; if ( !pamdbpsw ) pamdbpsw = ""; if ( strcmp(pamdbhost,"") ) host = pamdbhost; if ( strcmp(pamdbuser,"") ) user = pamdbuser; if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw; // // TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data()); if ( !dbc->IsConnected() ) return 1; stringstream myquery; myquery.str(""); myquery << "SET time_zone='+0:00'"; dbc->Query(myquery.str().c_str()); GL_PARAM *glparam = new GL_PARAM(); glparam->Query_GL_PARAM(1,1,dbc); // parameters stored in DB in GL_PRAM table trk->LoadField(glparam->PATH+glparam->NAME); // Bool_t defcal = true; Int_t error=glparam->Query_GL_PARAM(run->RUNHEADER_TIME,201,dbc); // parameters stored in DB in GL_PRAM table if ( error<0 ) { return(1); }; printf(" Reading ToF parameter file: %s \n",(glparam->PATH+glparam->NAME).Data()); if ( (UInt_t)glparam->TO_TIME != (UInt_t)4294967295UL ) defcal = false; // Int_t nlen = (Int_t)(glparam->PATH+glparam->NAME).Length(); rdtofcal((char *)(glparam->PATH+glparam->NAME).Data(),&nlen); // Int_t adc[4][12]; Int_t tdc[4][12]; Float_t tdcc[4][12]; // // process tof data // for (Int_t hh=0; hh<12;hh++){ for (Int_t kk=0; kk<4;kk++){ adc[kk][hh] = 4095; tdc[kk][hh] = 4095; tdcc[kk][hh] = 4095.; tofinput_.adc[hh][kk] = 4095; tofinput_.tdc[hh][kk] = 4095; }; }; Int_t ntrkentry = 0; Int_t npmtentry = 0; Int_t gg = 0; Int_t hh = 0; Int_t adcf[48]; memset(adcf, 0, 48*sizeof(Int_t)); Int_t tdcf[48]; memset(tdcf, 0, 48*sizeof(Int_t)); for (Int_t pm=0; pm < this->ntrk() ; pm++){ ToFTrkVar *ttf = this->GetToFTrkVar(pm); for ( Int_t nc=0; nc < ttf->npmttdc; nc++){ if ( (ttf->tdcflag).At(nc) != 0 ) tdcf[(ttf->pmttdc).At(nc)] = 1; }; for ( Int_t nc=0; nc < ttf->npmtadc; nc++){ if ( (ttf->adcflag).At(nc) != 0 ) adcf[(ttf->pmtadc).At(nc)] = 1; }; }; // for (Int_t pm=0; pm < this->npmt() ; pm++){ ToFPMT *pmt = this->GetToFPMT(pm); this->GetPMTIndex(pmt->pmt_id, gg, hh); if ( adcf[pmt->pmt_id] == 0 ){ tofinput_.adc[gg][hh] = (int)pmt->adc; adc[hh][gg] = (int)pmt->adc; }; if ( tdcf[pmt->pmt_id] == 0 ){ tofinput_.tdc[gg][hh] = (int)pmt->tdc; tdc[hh][gg] = (int)pmt->tdc; }; tdcc[hh][gg] = (float)pmt->tdc_tw; // Int_t pppid = this->GetPMTid(hh,gg); // 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); }; // Int_t unpackError = this->unpackError; // for (Int_t hh=0; hh<5;hh++){ tofinput_.patterntrig[hh]=trg->patterntrig[hh]; }; // this->Clear(); // Int_t pmt_id = 0; ToFPMT *t_pmt = new ToFPMT(); if(!(this->PMT)) this->PMT = new TClonesArray("ToFPMT",12); //ELENA TClonesArray &tpmt = *this->PMT; ToFTrkVar *t_tof = new ToFTrkVar(); if(!(this->ToFTrk)) this->ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA TClonesArray &t = *this->ToFTrk; // // // Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract common and track-related variables. // npmtentry = 0; // ntrkentry = 0; // // Calculate tracks informations from ToF alone // tofl2com(); // memcpy(this->tof_j_flag,tofoutput_.tof_j_flag,6*sizeof(Int_t)); // t_tof->trkseqno = -1; // // and now we must copy from the output structure to the level2 class: // t_tof->npmttdc = 0; // for (Int_t hh=0; hh<12;hh++){ for (Int_t kk=0; kk<4;kk++){ if ( tofoutput_.tofmask[hh][kk] != 0 ){ pmt_id = this->GetPMTid(kk,hh); t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc); t_tof->tdcflag.AddAt(tofoutput_.tdcflagtof[hh][kk],t_tof->npmttdc); // gf: Jan 09/07 t_tof->npmttdc++; }; }; }; for (Int_t kk=0; kk<13;kk++){ t_tof->beta[kk] = tofoutput_.betatof_a[kk]; } // t_tof->npmtadc = 0; for (Int_t hh=0; hh<12;hh++){ for (Int_t kk=0; kk<4;kk++){ if ( tofoutput_.adctof_c[hh][kk] < 1000 ){ t_tof->dedx.AddAt(tofoutput_.adctof_c[hh][kk],t_tof->npmtadc); pmt_id = this->GetPMTid(kk,hh); t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc); t_tof->adcflag.AddAt(tofoutput_.adcflagtof[hh][kk],t_tof->npmtadc); // gf: Jan 09/07 t_tof->npmtadc++; }; }; }; // memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos)); memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos)); memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof)); memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof)); // new(t[ntrkentry]) ToFTrkVar(*t_tof); ntrkentry++; t_tof->Clear(); // // // t_pmt->Clear(); // for (Int_t hh=0; hh<12;hh++){ for (Int_t kk=0; kk<4;kk++){ // new WM if ( tofoutput_.tdc_c[hh][kk] < 4095 || adc[kk][hh] < 4095 || tdc[kk][hh] < 4095 ){ // if ( tdcc[kk][hh] < 4095. || adc[kk][hh] < 4095 || tdc[kk][hh] < 4095 ){ // t_pmt->pmt_id = this->GetPMTid(kk,hh); t_pmt->tdc_tw = tofoutput_.tdc_c[hh][kk]; t_pmt->adc = (Float_t)adc[kk][hh]; t_pmt->tdc = (Float_t)tdc[kk][hh]; // new(tpmt[npmtentry]) ToFPMT(*t_pmt); npmtentry++; t_pmt->Clear(); }; }; }; // // Calculate track-related variables // if ( trk->ntrk() > 0 ){ // // We have at least one track // // // Run over tracks // for(Int_t nt=0; nt < trk->ntrk(); nt++){ // TrkTrack *ptt = trk->GetStoredTrack(nt); // // Copy the alpha vector in the input structure // for (Int_t e = 0; e < 5 ; e++){ tofinput_.al_pp[e] = ptt->al[e]; }; // // Get tracker related variables for this track // toftrk(); // // Copy values in the class from the structure (we need to use a temporary class to store variables). // t_tof->npmttdc = 0; for (Int_t hh=0; hh<12;hh++){ for (Int_t kk=0; kk<4;kk++){ if ( tofoutput_.tofmask[hh][kk] != 0 ){ pmt_id = this->GetPMTid(kk,hh); t_tof->pmttdc.AddAt(pmt_id,t_tof->npmttdc); t_tof->tdcflag.AddAt(tofoutput_.tdcflag[hh][kk],t_tof->npmttdc); // gf: Jan 09/07 t_tof->npmttdc++; }; }; }; for (Int_t kk=0; kk<13;kk++){ t_tof->beta[kk] = tofoutput_.beta_a[kk]; }; // t_tof->npmtadc = 0; for (Int_t hh=0; hh<12;hh++){ for (Int_t kk=0; kk<4;kk++){ if ( tofoutput_.adc_c[hh][kk] < 1000 ){ t_tof->dedx.AddAt(tofoutput_.adc_c[hh][kk],t_tof->npmtadc); pmt_id = this->GetPMTid(kk,hh); t_tof->pmtadc.AddAt(pmt_id,t_tof->npmtadc); t_tof->adcflag.AddAt(tofoutput_.adcflag[hh][kk],t_tof->npmtadc); // gf: Jan 09/07 t_tof->npmtadc++; }; }; }; // memcpy(t_tof->xtofpos,tofoutput_.xtofpos,sizeof(t_tof->xtofpos)); memcpy(t_tof->ytofpos,tofoutput_.ytofpos,sizeof(t_tof->ytofpos)); memcpy(t_tof->xtr_tof,tofoutput_.xtr_tof,sizeof(t_tof->xtr_tof)); memcpy(t_tof->ytr_tof,tofoutput_.ytr_tof,sizeof(t_tof->ytr_tof)); // // Store the tracker track number in order to be sure to have shyncronized data during analysis // t_tof->trkseqno = nt; // // create a new object for this event with track-related variables // new(t[ntrkentry]) ToFTrkVar(*t_tof); ntrkentry++; t_tof->Clear(); // }; // loop on all the tracks // this->unpackError = unpackError; if ( defcal ){ this->default_calib = 1; } else { this->default_calib = 0; }; }; return(0); } ToFdEdx::ToFdEdx() { memset(conn,0,12*sizeof(Bool_t)); memset(ts,0,12*sizeof(UInt_t)); memset(te,0,12*sizeof(UInt_t)); Define_PMTsat(); Clear(); } //------------------------------------------------------------------------ void ToFdEdx::CheckConnectors(UInt_t atime, GL_PARAM *glparam, TSQLServer *dbc) { for(int i=0; i<12; i++){ if(atime<=ts[i] || atime>te[i]){ Int_t error=glparam->Query_GL_PARAM(atime,210+i,dbc); // parameters stored in DB in GL_PRAM table if ( error<0 ) { conn[i]=false; ts[i]=0; te[i]=numeric_limits::max(); }; if ( !error ){ conn[i]=true; ts[i]=glparam->FROM_TIME; te[i]=glparam->TO_TIME; } if ( error>0 ){ conn[i]=false; ts[i]=glparam->TO_TIME; TSQLResult *pResult; TSQLRow *row; 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); pResult=dbc->Query(query.Data()); if(!pResult->GetRowCount()){ te[i]=numeric_limits::max(); }else{ row=pResult->Next(); te[i]=(UInt_t)atoll(row->GetField(0)); } } // } } } //------------------------------------------------------------------------ void ToFdEdx::Clear(Option_t *option) { // // Set arrays and initialize structure eDEDXpmt.Set(48); eDEDXpmt.Reset(-1); // Set array size and reset structure eZpmt.Set(48); eZpmt.Reset(-1); eDEDXpad.Set(24); eDEDXpad.Reset(-1); eZpad.Set(24); eZpad.Reset(-1); eDEDXlayer.Set(6); eDEDXlayer.Reset(-1); eZlayer.Set(6); eZlayer.Reset(-1); eDEDXplane.Set(3); eDEDXplane.Reset(-1); eZplane.Set(3); eZplane.Reset(-1); INFOpmt.Set(48); INFOpmt.Reset(0); INFOlayer.Set(6); INFOlayer.Reset(0); // }; //------------------------------------------------------------------------ void ToFdEdx::Print(Option_t *option) { // printf("========================================================================\n"); }; //------------------------------------------------------------------------ // void ToFdEdx::InitPar(TString parname, TString parfile) // { // // expensive function - call it once/run // ReadParAtt( Form("%s/attenuation.txt" , pardir) ); // ReadParPos( Form("%s/desaturation_position.txt" , pardir) ); // ReadParBBneg( Form("%s/BetheBloch.txt" , pardir) ); // ReadParBBpos( Form("%s/BetheBloch_betagt1.txt" , pardir) ); // ReadParDesatBB( Form("%s/desaturation_beta.txt" , pardir) ); // }; //------------------------------------------------------------------------ void ToFdEdx::Process(UInt_t atime, Float_t betamean, Float_t *xtr_tof, Float_t *ytr_tof, pamela::tof::TofEvent *tofl0 ) { // the parameters should be already initialised by InitPar() Clear(); // Float_t betamean = fabs(trackTRK->GetToFTrack()->beta[12]); if(betamean<0.05 || betamean>2){ for(int i=0;i<48;i++)INFOpmt[i]=1; } // define angle: double dx = xtr_tof[1] - xtr_tof[5]; double dy = ytr_tof[0] - ytr_tof[4]; double dr = sqrt(dx*dx+dy*dy); double theta=atan(dr/76.81); // TArrayF adc; Float_t adc[48]; ToFLevel2 tf; for (Int_t gg=0; gg<4;gg++){ for (Int_t hh=0; hh<12;hh++){ // tofinput_.tdc[hh][gg]=tofEvent->tdc[gg][hh]; int mm = tf.GetPMTid(gg,hh); adc[mm]=tofl0->adc[gg][hh]; }; }; for( int ii=0; ii<48; ii++ ) { if( adc[ii] >= PMTsat[ii]-5 ) continue; if( adc[ii] <= 0. ) continue; double adcpC = f_adcPC( adc[ii] ); // - adc conversion in pC double adccorr = adcpC*fabs(cos(theta)); if(adccorr<=0.) continue; //--------------------- TABLE OF PERIODS WITH HV PROBLEMS ---------------------------- int Aconn=conn[0]; // PMT 0,20,22,24 int Bconn=conn[1]; // PMT 6,12,26,34 int Cconn=conn[2]; // PMT 4,14,28,32 int Dconn=conn[3]; // PMT 2,8,10,30 int Econn=conn[4]; // PMT 42,43,44,47 int Fconn=conn[5]; // PMT 7,19,23,27 int Gconn=conn[6]; // PMT 3,11,25,33 int Hconn=conn[7]; // PMT 1,9,13,21 int Iconn=conn[8]; // PMT 5,29,31,35 int Lconn=conn[9]; // PMT 37,40,45,46 int Mconn=conn[10]; // PMT 15,16,17,18 int Nconn=conn[11]; // PMT 36,38,39,41 // int standard=0; if( false ) cout << Gconn << Iconn << Lconn <=1153660001 && atime<=1154375000)Dconn=1; // else if(atime>=1155850001 && atime<=1156280000){ // Hconn=1; // Nconn=1; // } // else if(atime>=1168490001 && atime<=1168940000)Dconn=1; // else if(atime>=1168940001 && atime<=1169580000){ // Fconn=1; // Mconn=1; // } // else if(atime>=1174665001 && atime<=1175000000)Bconn=1; // else if(atime>=1176120001 && atime<=1176800000)Hconn=1; // else if(atime>=1176800001 && atime<=1178330000)Econn=1; // else if(atime>=1178330001 && atime<=1181322000)Hconn=1; // else if(atime>=1182100001 && atime<=1183030000)Aconn=1; // else if(atime>=1184000001 && atime<=1184570000)Hconn=1; // else if(atime>=1185090001 && atime<=1185212000)Dconn=1; // else if(atime>=1191100001 && atime<=1191940000)Dconn=1; // else if(atime>=1196230001 && atime<=1196280000)Hconn=1; // else if(atime>=1206100001 && atime<=1206375600)Cconn=1; // else if(atime>=1217989201 && atime<=1218547800)Econn=1; // else if(atime>=1225789201 && atime<=1226566800)Econn=1; // else if(atime>=1229400901 && atime<=1229700000)Econn=1; // else if(atime>=1230318001 && atime<=1230415200)Econn=1; // else { // standard=1; // } if(atime<1158720000)S115B_ok=1; else S115B_break=1; //------------------------------------------------------------------------ //---------------------------------------------------- Z reconstruction double adcHe, adcnorm, adclin, dEdx, Zeta; adcHe=-2; adcnorm=-2; adclin=-2; dEdx=-2; Zeta=-2; // float ZetaH=-2; // float dEdxH=-2; // double day = (atime-1150000000)/84600; if(Aconn==1 && (ii==0 || ii==20 || ii==22 || ii==24)){ adcHe = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.675; } else if(Bconn==1 && (ii==6 || ii==12 || ii==26 || ii==34)){ adcHe = (Get_adc_he(ii, xtr_tof, ytr_tof))/2.482; } else if(Cconn==1 && (ii==4 || ii==14 || ii==28 || ii==32)){ adcHe = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.464; } else if(Dconn==1 && (ii==2 || ii==8 || ii==10 || ii==30)){ adcHe = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.995; } else if(Econn==1 && (ii==42 || ii==43 || ii==44 || ii==47)){ adcHe = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.273; } else if(Fconn==1 && (ii==7 || ii==19 || ii==23 || ii==27)){ adcHe = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.565; } else if(Mconn==1 && (ii==15 || ii==16 || ii==17 || ii==18)){ adcHe = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.565; } else if(Nconn==1 && (ii==36 || ii==38 || ii==39 || ii==41)){ adcHe = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.018; } else if(Hconn==1 && (ii==1 || ii==13 || ii==21 || (ii==9&&S115B_ok==1))){ adcHe = (Get_adc_he(ii, xtr_tof, ytr_tof))/1.84; } else if(S115B_break==1 && ii==9 && Hconn==0){ adcHe = f_att5B( ytr_tof[0] ); //N.B.: this function refers to the Carbon!!! } else if(S115B_break==1 && ii==9 && Hconn==1){ adcHe = (f_att5B( ytr_tof[0] ))/1.64; } else adcHe = Get_adc_he(ii, xtr_tof, ytr_tof); if(adcHe<=0) continue; if(ii==9 && S115B_break==1) adcnorm = f_pos5B(adccorr); else adcnorm = f_pos( (parPos[ii]), adccorr); if(adcnorm<=0) continue; if(ii==9 && S115B_break==1) adclin = 36.*adcnorm/adcHe; else adclin = 4.*adcnorm/adcHe; if(adclin<=0) continue; double dEdxHe=-2; if(ii==9 && S115B_break==1){ if( betamean <1. ) dEdxHe = f_BB5B( betamean ); else dEdxHe = 33; } else { if( betamean <1. ) dEdxHe = f_BB( (parBBneg[ii]), betamean ); else dEdxHe = parBBpos[ii]; } if(dEdxHe<=0) continue; if(ii==9 && S115B_break==1) dEdx = f_desatBB5B( adclin ); else dEdx = f_desatBB((parDesatBB[ii]), adclin ); if(dEdx<=0) continue; if(ii==9 && S115B_break==1) Zeta = sqrt(36.*(dEdx/dEdxHe)); else Zeta = sqrt(4.*(dEdx/dEdxHe)); if(Zeta<=0) continue; //--------------------- TIME DEPENDENCE ---------------------------------- // TArrayF &binx = TDx[ii]; // TArrayF &biny = TDy[ii]; // for(int k=0; k<200; k++) { // if (day > binx[k]-2.5 && day<=binx[k]+2.5 && biny[k]>0) { // ZetaH=Zeta/biny[k]*6; // dEdxH=dEdx/(pow(biny[k],2))*36; // } // } // if(ZetaH!=-2)eZpmt[ii]=(Float_t)ZetaH; // else eZpmt[ii]=(Float_t)Zeta; // if(dEdxH!=-2)eDEDXpmt[ii]=(Float_t)dEdxH; // else eDEDXpmt[ii]=(Float_t)dEdx; // printf("%5d %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %5.4f \n", ii, adcpC, adccorr, adcHe, dEdxHe, dEdx, Zeta, betamean ); eZpmt[ii]=(Float_t)Zeta; eDEDXpmt[ii]=(Float_t)dEdx; } //end loop on 48 PMT //--------------------------------------------------- paddle + layer -------------------- for(int j=0;j<48;j++){ int k=100; if(j%2==0 || j==0)k=j/2; double zpdl=-1; if((j%2==0 || j==0) && eZpmt[j]!=-1 && eZpmt[j+1]!=-1){ zpdl=0.5*(eZpmt[j]+eZpmt[j+1]); }else if((j%2==0 || j==0) && eZpmt[j]!=-1 && eZpmt[j+1]==-1){ zpdl=eZpmt[j]; }else if((j%2==0 || j==0) && eZpmt[j]==-1 && eZpmt[j+1]!=-1){ zpdl=eZpmt[j+1]; } if(j%2==0 || j==0)eZpad[k]= (Float_t)zpdl; if((j%2==0 || j==0)&&eZpad[k]!=-1){ if(k>=0&&k<8)eZlayer[0]=eZpad[k]; if(k>=8&&k<14)eZlayer[1]=eZpad[k]; if(k>=14&&k<16)eZlayer[2]=eZpad[k]; if(k>=16&&k<18)eZlayer[3]=eZpad[k]; if(k>=18&&k<21)eZlayer[4]=eZpad[k]; if(k>=21)eZlayer[5]=eZpad[k]; } if(eZlayer[0]!=-1&&eZlayer[1]!=-1&&fabs(eZlayer[0]-eZlayer[1])<1.5)eZplane[0]=0.5*(eZlayer[0]+eZlayer[1]); else if(eZlayer[0]!=-1&&eZlayer[1]==-1)eZplane[0]=eZlayer[0]; else if(eZlayer[1]!=-1&&eZlayer[0]==-1)eZplane[0]=eZlayer[1]; if(eZlayer[2]!=-1&&eZlayer[3]!=-1&&fabs(eZlayer[2]-eZlayer[3])<1.5)eZplane[1]=0.5*(eZlayer[2]+eZlayer[3]); else if(eZlayer[2]!=-1&&eZlayer[3]==-1)eZplane[1]=eZlayer[2]; else if(eZlayer[3]!=-1&&eZlayer[2]==-1)eZplane[1]=eZlayer[3]; if(eZlayer[4]!=-1&&eZlayer[5]!=-1&&fabs(eZlayer[4]-eZlayer[5])<1.5)eZplane[2]=0.5*(eZlayer[4]+eZlayer[5]); else if(eZlayer[4]!=-1&&eZlayer[5]==-1)eZplane[2]=eZlayer[4]; else if(eZlayer[5]!=-1&&eZlayer[4]==-1)eZplane[2]=eZlayer[5]; } for(int jj=0;jj<48;jj++){ int k=100; if(jj%2==0 || jj==0)k=jj/2; double dedxpdl=-1; if((jj%2==0 || jj==0) && eDEDXpmt[jj]!=-1 && eDEDXpmt[jj+1]!=-1){ dedxpdl=0.5*(eDEDXpmt[jj]+eDEDXpmt[jj+1]); }else if((jj%2==0 || jj==0) && eDEDXpmt[jj]!=-1 && eDEDXpmt[jj+1]==-1){ dedxpdl=eDEDXpmt[jj]; }else if((jj%2==0 || jj==0) && eDEDXpmt[jj]==-1 && eDEDXpmt[jj+1]!=-1){ dedxpdl=eDEDXpmt[jj+1]; } if(jj%2==0 || jj==0)eDEDXpad[k]= (Float_t)dedxpdl; if((jj%2==0 || jj==0)&&eDEDXpad[k]!=-1){ if(k>=0&&k<8)eDEDXlayer[0]=eDEDXpad[k]; if(k>=8&&k<14)eDEDXlayer[1]=eDEDXpad[k]; if(k>=14&&k<16)eDEDXlayer[2]=eDEDXpad[k]; if(k>=16&&k<18)eDEDXlayer[3]=eDEDXpad[k]; if(k>=18&&k<21)eDEDXlayer[4]=eDEDXpad[k]; if(k>=21)eDEDXlayer[5]=eDEDXpad[k]; } if(eDEDXlayer[0]!=-1&&eDEDXlayer[1]!=-1&&fabs(eDEDXlayer[0]-eDEDXlayer[1])<10)eDEDXplane[0]=0.5*(eDEDXlayer[0]+eDEDXlayer[1]); else if(eDEDXlayer[0]!=-1&&eDEDXlayer[1]==-1)eDEDXplane[0]=eDEDXlayer[0]; else if(eDEDXlayer[1]!=-1&&eDEDXlayer[0]==-1)eDEDXplane[0]=eDEDXlayer[1]; if(eDEDXlayer[2]!=-1&&eDEDXlayer[3]!=-1&&fabs(eDEDXlayer[2]-eDEDXlayer[3])<10)eDEDXplane[1]=0.5*(eDEDXlayer[2]+eDEDXlayer[3]); else if(eDEDXlayer[2]!=-1&&eDEDXlayer[3]==-1)eDEDXplane[1]=eDEDXlayer[2]; else if(eDEDXlayer[3]!=-1&&eDEDXlayer[2]==-1)eDEDXplane[1]=eDEDXlayer[3]; if(eDEDXlayer[4]!=-1&&eDEDXlayer[5]!=-1&&fabs(eDEDXlayer[4]-eDEDXlayer[5])<10)eDEDXplane[2]=0.5*(eDEDXlayer[4]+eDEDXlayer[5]); else if(eDEDXlayer[4]!=-1&&eDEDXlayer[5]==-1)eDEDXplane[2]=eDEDXlayer[4]; else if(eDEDXlayer[5]!=-1&&eDEDXlayer[4]==-1)eDEDXplane[2]=eDEDXlayer[5]; } }; //------------------------------------------------------------------------ void ToFdEdx::PrintTD() { for(int i=0; i<48; i++) { TArrayF &binx = TDx[i]; TArrayF &biny = TDy[i]; for(int k=0; k<200; k++) { // bin temporali printf("%d %d %f %f", i,k, binx[k], biny[k]); } } } //------------------------------------------------------------------------ void ToFdEdx::Define_PMTsat() { Float_t sat[48] = { 3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27, 3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14, 3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62, 3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26, 3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18, 3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 }; PMTsat.Set(48,sat); } //------------------------------------------------------------------------ // void ToFdEdx::ReadParTD( Int_t ipmt, const char *fname ) // { // printf("read %s\n",fname); // if(ipmt<0) return; // if(ipmt>47) return; // FILE *fattin = fopen( fname , "r" ); // Float_t yTD[200],xTD[200]; // for(int j=0;j<200;j++){ // float x,y,ym,e; // if(fscanf(fattin,"%f %f %f %f", // &x, &y, &ym, &e )!=4) break; // xTD[j]=x; // if(ym>0&&fabs(y-ym)>1) yTD[j]=ym; // else yTD[j]=y; // } // TDx[ipmt].Set(200,xTD); // TDy[ipmt].Set(200,yTD); // fclose(fattin); // } //------------------------------------------------------------------------ void ToFdEdx::ReadParBBpos( const char *fname ) { printf("read %s\n",fname); parBBpos.Set(48); FILE *fattin = fopen( fname , "r" ); for (int i=0; i<48; i++) { int tid=0; float tp; if(fscanf(fattin,"%d %f", &tid, &tp )!=2) break; parBBpos[i]=tp; } fclose(fattin); } //------------------------------------------------------------------------ void ToFdEdx::ReadParDesatBB( const char *fname ) { printf("read %s\n",fname); FILE *fattin = fopen( fname , "r" ); for (int i=0; i<48; i++) { int tid=0; float tp[3]; if(fscanf(fattin,"%d %f %f %f", &tid, &tp[0], &tp[1], &tp[2] )!=4) break; parDesatBB[i].Set(3,tp); } fclose(fattin); } //------------------------------------------------------------------------ void ToFdEdx::ReadParBBneg( const char *fname ) { printf("read %s\n",fname); FILE *fattin = fopen( fname , "r" ); for (int i=0; i<48; i++) { int tid=0; float tp[3]; if(fscanf(fattin,"%d %f %f %f", &tid, &tp[0], &tp[1], &tp[2] )!=4) break; parBBneg[i].Set(3,tp); } fclose(fattin); } //------------------------------------------------------------------------ void ToFdEdx::ReadParPos( const char *fname ) { printf("read %s\n",fname); FILE *fattin = fopen( fname , "r" ); for (int i=0; i<48; i++) { int tid=0; float tp[4]; if(fscanf(fattin,"%d %f %f %f %f", &tid, &tp[0], &tp[1], &tp[2], &tp[3])!=5) break; parPos[i].Set(4,tp); } fclose(fattin); } //------------------------------------------------------------------------ void ToFdEdx::ReadParAtt( const char *fname ) { printf("read %s\n",fname); FILE *fattin = fopen( fname , "r" ); for (int i=0; i<48; i++) { int tid=0; float tp[6]; if(fscanf(fattin,"%d %f %f %f %f %f %f", &tid, &tp[0], &tp[1], &tp[2], &tp[3], &tp[4], &tp[5] )!=7) break; parAtt[i].Set(6,tp); } fclose(fattin); } double ToFdEdx::f_att( TArrayF &p, float x ) { return p[0] + p[1]*x + p[2]*x*x + p[3]*x*x*x + p[4]*x*x*x*x + p[5]*x*x*x*x*x; } //------------------------------------------------------------------------ double ToFdEdx::f_att5B( float x ) { return 101.9409 + 6.643781*x + 0.2765518*x*x + 0.004617647*x*x*x + 0.0006195132*x*x*x*x + 0.00002813734*x*x*x*x*x; } double ToFdEdx::f_pos( TArrayF &p, float x ) { return p[0] + p[1]*x + p[2]*x*x + p[3]*x*x*x; } double ToFdEdx::f_pos5B( float x ) { return 15.45132 + 0.8369721*x + 0.0005*x*x; } double ToFdEdx::f_adcPC( float x ) { return 28.12+0.6312*x-5.647e-05*x*x+3.064e-08*x*x*x; } float ToFdEdx::Get_adc_he( int id, float pl_x[6], float pl_y[6]) { // // input: id - pmt [0:47} // pl_x - coord x of the tof plane // pl_y - coord y adc_he = 0; if( eGeom.GetXY(id)==1 ) adc_he = f_att( (parAtt[id]), pl_x[eGeom.GetPlane(id)] ); if( eGeom.GetXY(id)==2 ) adc_he = f_att( (parAtt[id]), pl_y[eGeom.GetPlane(id)] ); return adc_he; } //------------------------------------------------------------------------ double ToFdEdx::f_BB( TArrayF &p, float x ) { return p[0]/(x*x)*(log(x*x/(1-x*x)) - p[1]*x*x - p[2]); } //------------------------------------------------------------------------ double ToFdEdx::f_BB5B( float x ) { return 0.165797/(x*x)*(log(x*x/(1-x*x)) + 140.481*x*x + 52.9258); } //------------------------------------------------------------------------ double ToFdEdx::f_desatBB( TArrayF &p, float x ) { return p[0] + p[1]*x + p[2]*x*x; } //------------------------------------------------------------------------ double ToFdEdx::f_desatBB5B( float x ) { return -2.4 + 0.75*x + 0.009*x*x; }