diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx index 58c076f5954521d3c75520caaa84e4df4b1130ba..6eb6bbd2ed91291fa01f7fbd600b1e187dfafcd0 100644 --- a/src/THcHodoscope.cxx +++ b/src/THcHodoscope.cxx @@ -293,6 +293,17 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date ) fHodoPosInvAdcAdc=new Double_t [fMaxHodoScin]; fHodoNegInvAdcAdc=new Double_t [fMaxHodoScin]; + //New Time-Walk Calibration Parameters + fHodoVelFit=new Double_t [fMaxHodoScin]; + fHodoCableFit=new Double_t [fMaxHodoScin]; + fHodo_LCoeff=new Double_t [fMaxHodoScin]; + fHodoPos_c1=new Double_t [fMaxHodoScin]; + fHodoNeg_c1=new Double_t [fMaxHodoScin]; + fHodoPos_c2=new Double_t [fMaxHodoScin]; + fHodoNeg_c2=new Double_t [fMaxHodoScin]; + fHodoSigmaPos=new Double_t [fMaxHodoScin]; + fHodoSigmaNeg=new Double_t [fMaxHodoScin]; + fNHodoscopes = 2; fxLoScin = new Int_t [fNHodoscopes]; fxHiScin = new Int_t [fNHodoscopes]; @@ -322,10 +333,10 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date ) {"scin_tdc_max", &fScinTdcMax, kDouble}, {"tof_tolerance", &fTofTolerance, kDouble, 0, 1}, {"pathlength_central", &fPathLengthCentral, kDouble}, - {"hodo_pos_sigma", &fHodoPosSigma[0], kDouble, fMaxHodoScin}, - {"hodo_neg_sigma", &fHodoNegSigma[0], kDouble, fMaxHodoScin}, - {"hodo_pos_ped_limit", &fHodoPosPedLimit[0], kInt, fMaxHodoScin}, - {"hodo_neg_ped_limit", &fHodoNegPedLimit[0], kInt, fMaxHodoScin}, + {"hodo_pos_sigma", &fHodoPosSigma[0], kDouble, fMaxHodoScin, 1}, + {"hodo_neg_sigma", &fHodoNegSigma[0], kDouble, fMaxHodoScin, 1}, + {"hodo_pos_ped_limit", &fHodoPosPedLimit[0], kInt, fMaxHodoScin, 1}, + {"hodo_neg_ped_limit", &fHodoNegPedLimit[0], kInt, fMaxHodoScin, 1}, {"tofusinginvadc", &fTofUsingInvAdc, kInt, 0, 1}, {"xloscin", &fxLoScin[0], kInt, (UInt_t) fNHodoscopes}, {"xhiscin", &fxHiScin[0], kInt, (UInt_t) fNHodoscopes}, @@ -359,6 +370,11 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date ) fHodoPosAdcTimeWindowMax[ip] = 1000.; fHodoNegAdcTimeWindowMin[ip] = -1000.; fHodoNegAdcTimeWindowMax[ip] = 1000.; + + fHodoPosPedLimit[ip] = 0.0; + fHodoNegPedLimit[ip] = 0.0; + fHodoPosSigma[ip] = 0.2; + fHodoNegSigma[ip] = 0.2; } fTOFCalib_shtrk_lo=-kBig; fTOFCalib_shtrk_hi= kBig; @@ -410,7 +426,7 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date ) if (fTofUsingInvAdc) { DBRequest list2[]={ - {"hodo_vel_light", &fHodoVelLight[0], kDouble, fMaxHodoScin}, + {"hodo_vel_light", &fHodoVelLight[0], kDouble, fMaxHodoScin, 1}, {"hodo_pos_invadc_offset",&fHodoPosInvAdcOffset[0],kDouble,fMaxHodoScin}, {"hodo_neg_invadc_offset",&fHodoNegInvAdcOffset[0],kDouble,fMaxHodoScin}, {"hodo_pos_invadc_linear",&fHodoPosInvAdcLinear[0],kDouble,fMaxHodoScin}, @@ -419,10 +435,18 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date ) {"hodo_neg_invadc_adc",&fHodoNegInvAdcAdc[0],kDouble,fMaxHodoScin}, {0} }; + + for (UInt_t i=0; i<fMaxHodoScin; i++) + { + //Set scin Velocity/Cable to default + fHodoVelLight[i] = 15.0; + + } + gHcParms->LoadParmValues((DBRequest*)&list2,prefix); }; - if (!fTofUsingInvAdc) { - DBRequest list3[]={ + /* if (!fTofUsingInvAdc) { + DBRequest list3[]={ {"hodo_vel_light", &fHodoVelLight[0], kDouble, fMaxHodoScin}, {"hodo_pos_minph", &fHodoPosMinPh[0], kDouble, fMaxHodoScin}, {"hodo_neg_minph", &fHodoNegMinPh[0], kDouble, fMaxHodoScin}, @@ -430,10 +454,49 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date ) {"hodo_neg_phc_coeff", &fHodoNegPhcCoeff[0], kDouble, fMaxHodoScin}, {"hodo_pos_time_offset", &fHodoPosTimeOffset[0], kDouble, fMaxHodoScin}, {"hodo_neg_time_offset", &fHodoNegTimeOffset[0], kDouble, fMaxHodoScin}, - {0} + {0} }; + + gHcParms->LoadParmValues((DBRequest*)&list3,prefix); - }; + */ + DBRequest list4[]={ + {"hodo_velFit", &fHodoVelFit[0], kDouble, fMaxHodoScin, 1}, + {"hodo_cableFit", &fHodoCableFit[0], kDouble, fMaxHodoScin, 1}, + {"hodo_LCoeff", &fHodo_LCoeff[0], kDouble, fMaxHodoScin, 1}, + {"c1_Pos", &fHodoPos_c1[0], kDouble, fMaxHodoScin, 1}, + {"c1_Neg", &fHodoNeg_c1[0], kDouble, fMaxHodoScin, 1}, + {"c2_Pos", &fHodoPos_c2[0], kDouble, fMaxHodoScin, 1}, + {"c2_Neg", &fHodoNeg_c2[0], kDouble, fMaxHodoScin, 1}, + {"TDC_threshold", &fTdc_Thrs, kDouble, 0, 1}, + {"hodo_PosSigma", &fHodoSigmaPos[0], kDouble, fMaxHodoScin, 1}, + {"hodo_NegSigma", &fHodoSigmaNeg[0], kDouble, fMaxHodoScin, 1}, + {0} + }; + + fTdc_Thrs = 1.0; + //Set Default Values if NOT defined in param file + for (UInt_t i=0; i<fMaxHodoScin; i++) + { + + //Turn OFF Time-Walk Correction if param file NOT found + fHodoPos_c1[i] = 0.0; + fHodoPos_c2[i] = 0.0; + fHodoNeg_c1[i] = 0.0; + fHodoNeg_c2[i] = 0.0; + } + for (UInt_t i=0; i<fMaxHodoScin; i++) + { + //Set scin Velocity/Cable to default + fHodoCableFit[i] = 0.0; + fHodoVelFit[i] = 15.0; + //set time coeff between paddles to default + fHodo_LCoeff[i] = 0.0; + + } + + gHcParms->LoadParmValues((DBRequest*)&list4,prefix); + if (fDebug >=1) { cout <<"******* Testing Hodoscope Parameter Reading ***\n"; cout<<"StarTimeCenter = "<<fStartTimeCenter<<endl; @@ -555,7 +618,16 @@ void THcHodoscope::DeleteArrays() delete [] fHodoNegAdcTimeWindowMax; fHodoNegAdcTimeWindowMax = NULL; delete [] fHodoPosAdcTimeWindowMin; fHodoPosAdcTimeWindowMin = NULL; delete [] fHodoPosAdcTimeWindowMax; fHodoPosAdcTimeWindowMax = NULL; - + + delete [] fHodoVelFit; fHodoVelFit = NULL; + delete [] fHodoCableFit; fHodoCableFit = NULL; + delete [] fHodo_LCoeff; fHodo_LCoeff = NULL; + delete [] fHodoPos_c1; fHodoPos_c1 = NULL; + delete [] fHodoNeg_c1; fHodoNeg_c1 = NULL; + delete [] fHodoPos_c2; fHodoPos_c2 = NULL; + delete [] fHodoNeg_c2; fHodoNeg_c2 = NULL; + delete [] fHodoSigmaPos; fHodoSigmaPos = NULL; + delete [] fHodoSigmaNeg; fHodoSigmaNeg = NULL; } //_____________________________________________________________________________ @@ -787,12 +859,21 @@ void THcHodoscope::EstimateFocalPlaneTime() if(twogoodtimes[ihhit]){ - Double_t sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) + - TMath::Power( fHodoNegSigma[GetScinIndex(ip,index)],2) ) ); + Double_t sigma = 0.0; + if(fTofUsingInvAdc) + { + sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) + + TMath::Power( fHodoNegSigma[GetScinIndex(ip,index)],2) ) ); + } + else{ + sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoSigmaPos[GetScinIndex(ip,index)],2) + + TMath::Power( fHodoSigmaNeg[GetScinIndex(ip,index)],2) ) ); + } + Double_t scinWeight = 1 / TMath::Power(sigma,2); Double_t zPosition = fPlanes[ip]->GetZpos() + (index%2)*fPlanes[ip]->GetDzpos(); // cout << "hit = " << ihhit + 1 << " zpos = " << zPosition << " sigma = " << sigma << endl; - + //cout << "fHodoSigma+ = " << fHodoSigmaPos[GetScinIndex(ip,index)] << endl; sumW += scinWeight; sumT += scinWeight * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime(); sumZ += scinWeight * zPosition; @@ -825,8 +906,17 @@ void THcHodoscope::EstimateFocalPlaneTime() Double_t zPosition = fPlanes[ip]->GetZpos() + (index%2)*fPlanes[ip]->GetDzpos(); Double_t timeDif = ( ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime() - t0 ); - Double_t sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) + - TMath::Power( fHodoNegSigma[GetScinIndex(ip,index)],2) ) ); + + Double_t sigma = 0.0; + if(fTofUsingInvAdc){ + sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) + + TMath::Power( fHodoNegSigma[GetScinIndex(ip,index)],2) ) ); + } + else { + sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoSigmaPos[GetScinIndex(ip,index)],2) + + TMath::Power( fHodoSigmaNeg[GetScinIndex(ip,index)],2) ) ); + } + fBetaNoTrkChiSq += ( ( zPosition / fBetaNoTrk - timeDif ) * ( zPosition / fBetaNoTrk - timeDif ) ) / ( sigma * sigma ); @@ -873,6 +963,7 @@ Double_t THcHodoscope::TimeWalkCorrection(const Int_t& paddle, Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) { + Int_t ntracks = tracks.GetLast()+1; // Number of reconstructed tracks // ------------------------------------------------- @@ -1000,6 +1091,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) Double_t tdc_pos = hit->GetPosTDC(); if(tdc_pos >=fScinTdcMin && tdc_pos <= fScinTdcMax ) { Double_t adc_pos = hit->GetPosADC(); + Double_t adcamp_pos = hit->GetPosADCpeak(); Double_t pathp = fPlanes[ip]->GetPosLeft() - scinLongCoord; fTOFPInfo[ihhit].pathp = pathp; Double_t timep = tdc_pos*fScinTdcToTime; @@ -1009,10 +1101,9 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) + fHodoPosInvAdcAdc[fPIndex] /TMath::Sqrt(TMath::Max(20.0*.020,adc_pos)); } else { - timep -= fHodoPosPhcCoeff[fPIndex]* - TMath::Sqrt(TMath::Max(0.0,adc_pos/fHodoPosMinPh[fPIndex]-1.0)) - + pathp/fHodoVelLight[fPIndex] - + fHodoPosTimeOffset[fPIndex]; + //Double_t tw_corr_pos = fHodoPos_c1[fPIndex]/pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) - fHodoPos_c1[fPIndex]/pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]); + Double_t tw_corr_pos = 1./pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) - 1./pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]); + timep += -tw_corr_pos + fHodo_LCoeff[fPIndex]; } fTOFPInfo[ihhit].scin_pos_time = timep; timep -= zcor; @@ -1023,6 +1114,7 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) Double_t tdc_neg = hit->GetNegTDC(); if(tdc_neg >=fScinTdcMin && tdc_neg <= fScinTdcMax ) { Double_t adc_neg = hit->GetNegADC(); + Double_t adcamp_neg = hit->GetNegADCpeak(); Double_t pathn = scinLongCoord - fPlanes[ip]->GetPosRight(); fTOFPInfo[ihhit].pathn = pathn; Double_t timen = tdc_neg*fScinTdcToTime; @@ -1032,10 +1124,10 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) + fHodoNegInvAdcAdc[fPIndex] /TMath::Sqrt(TMath::Max(20.0*.020,adc_neg)); } else { - timen -= fHodoNegPhcCoeff[fPIndex]* - TMath::Sqrt(TMath::Max(0.0,adc_neg/fHodoNegMinPh[fPIndex]-1.0)) - + pathn/fHodoVelLight[fPIndex] - + fHodoNegTimeOffset[fPIndex]; + // Double_t tw_corr_neg = fHodoNeg_c1[fPIndex]/pow(adcamp_neg/fTdc_Thrs,fHodoNeg_c2[fPIndex]) - fHodoNeg_c1[fPIndex]/pow(200./fTdc_Thrs, fHodoNeg_c2[fPIndex]); + Double_t tw_corr_neg = 1./pow(adcamp_neg/fTdc_Thrs,fHodoNeg_c2[fPIndex]) - 1./pow(200./fTdc_Thrs, fHodoNeg_c2[fPIndex]); + timen += -tw_corr_neg- 2*fHodoCableFit[fPIndex] + fHodo_LCoeff[fPIndex]; + } fTOFPInfo[ihhit].scin_neg_time = timen; timen -= zcor; @@ -1118,14 +1210,29 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) fTOFPInfo[ih].scin_neg_time ) / 2.; fTOFCalc[ih].scin_time_fp = ( fTOFPInfo[ih].time_pos + fTOFPInfo[ih].time_neg ) / 2.; - fTOFCalc[ih].scin_sigma = TMath::Sqrt( fHodoPosSigma[fPIndex] * fHodoPosSigma[fPIndex] + - fHodoNegSigma[fPIndex] * fHodoNegSigma[fPIndex] )/2.; + + if (fTofUsingInvAdc){ + fTOFCalc[ih].scin_sigma = TMath::Sqrt( fHodoPosSigma[fPIndex] * fHodoPosSigma[fPIndex] + + fHodoNegSigma[fPIndex] * fHodoNegSigma[fPIndex] )/2.; + } + else { + fTOFCalc[ih].scin_sigma = TMath::Sqrt( fHodoSigmaPos[fPIndex] * fHodoSigmaPos[fPIndex] + + fHodoSigmaNeg[fPIndex] * fHodoSigmaNeg[fPIndex] )/2.; + } + fTOFCalc[ih].good_scin_time = kTRUE; fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE; } else{ fTOFCalc[ih].scin_time = fTOFPInfo[ih].scin_pos_time; fTOFCalc[ih].scin_time_fp = fTOFPInfo[ih].time_pos; - fTOFCalc[ih].scin_sigma = fHodoPosSigma[fPIndex]; + + if (fTofUsingInvAdc){ + fTOFCalc[ih].scin_sigma = fHodoPosSigma[fPIndex]; + } + else{ + fTOFCalc[ih].scin_sigma = fHodoSigmaPos[fPIndex]; + } + fTOFCalc[ih].good_scin_time = kTRUE; fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE; } @@ -1133,7 +1240,12 @@ Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks ) if ( fTOFCalc[ih].good_tdc_neg ){ fTOFCalc[ih].scin_time = fTOFPInfo[ih].scin_neg_time; fTOFCalc[ih].scin_time_fp = fTOFPInfo[ih].time_neg; - fTOFCalc[ih].scin_sigma = fHodoNegSigma[fPIndex]; + if (fTofUsingInvAdc){ + fTOFCalc[ih].scin_sigma = fHodoNegSigma[fPIndex]; + } + else{ + fTOFCalc[ih].scin_sigma = fHodoSigmaNeg[fPIndex]; + } fTOFCalc[ih].good_scin_time = kTRUE; fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE; } diff --git a/src/THcHodoscope.h b/src/THcHodoscope.h index c96e0ede9cea7b74fb0ec29d3b2e8e79e97f0504..d9d83525e21dd3512744859b2fbc3ef38b5f6185 100644 --- a/src/THcHodoscope.h +++ b/src/THcHodoscope.h @@ -81,6 +81,18 @@ public: Double_t GetHodoPosAdcTimeWindowMin(Int_t iii) const {return fHodoPosAdcTimeWindowMin[iii];} Double_t GetHodoNegAdcTimeWindowMax(Int_t iii) const {return fHodoNegAdcTimeWindowMax[iii];} Double_t GetHodoNegAdcTimeWindowMin(Int_t iii) const {return fHodoNegAdcTimeWindowMin[iii];} + + //Get Time Walk Parameters + Double_t GetHodoVelFit(Int_t iii) const {return fHodoVelFit[iii];} + Double_t GetHodoCableFit(Int_t iii) const {return fHodoCableFit[iii];} + Double_t GetHodoLCoeff(Int_t iii) const {return fHodo_LCoeff[iii];} + + + Double_t GetHodoPos_c1(Int_t iii) const {return fHodoPos_c1[iii];} + Double_t GetHodoNeg_c1(Int_t iii) const {return fHodoNeg_c1[iii];} + Double_t GetHodoPos_c2(Int_t iii) const {return fHodoPos_c2[iii];} + Double_t GetHodoNeg_c2(Int_t iii) const {return fHodoNeg_c2[iii];} + Double_t GetTDCThrs() const {return fTdc_Thrs;} Double_t GetStartTimeCenter() const {return fStartTimeCenter;} Double_t GetStartTimeSlop() const {return fStartTimeSlop;} @@ -187,6 +199,18 @@ protected: Double_t* fHodoPosInvAdcAdc; Double_t* fHodoNegInvAdcAdc; + //New Time-Walk Calibration Parameters + Double_t* fHodoVelFit; + Double_t* fHodoCableFit; + Double_t* fHodo_LCoeff; + Double_t* fHodoPos_c1; + Double_t* fHodoNeg_c1; + Double_t* fHodoPos_c2; + Double_t* fHodoNeg_c2; + Double_t fTdc_Thrs; + Double_t* fHodoSigmaPos; + Double_t* fHodoSigmaNeg; + Double_t fPartMass; // Nominal particle mass Double_t fBetaNominal; // Beta for central ray of nominal particle type diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx index ef8f19a20ff820f2e5d582ff64b931c2a3e5958c..4f3896f771b0bc78877fc8dcdf68b205521c14e8 100644 --- a/src/THcScintillatorPlane.cxx +++ b/src/THcScintillatorPlane.cxx @@ -141,6 +141,14 @@ THcScintillatorPlane::~THcScintillatorPlane() delete [] fHodoNegInvAdcLinear; fHodoNegInvAdcLinear = NULL; delete [] fHodoPosInvAdcAdc; fHodoPosInvAdcAdc = NULL; delete [] fHodoNegInvAdcAdc; fHodoNegInvAdcAdc = NULL; + delete [] fHodoVelFit; fHodoVelFit = NULL; + delete [] fHodoCableFit; fHodoCableFit = NULL; + delete [] fHodo_LCoeff; fHodo_LCoeff = NULL; + delete [] fHodoPos_c1; fHodoPos_c1 = NULL; + delete [] fHodoNeg_c1; fHodoNeg_c1 = NULL; + delete [] fHodoPos_c2; fHodoPos_c2 = NULL; + delete [] fHodoNeg_c2; fHodoNeg_c2 = NULL; + delete [] fHodoVelLight; fHodoVelLight = NULL; delete [] fHodoSigma; fHodoSigma = NULL; @@ -270,6 +278,16 @@ Int_t THcScintillatorPlane::ReadDatabase( const TDatime& date ) fHodoPosInvAdcAdc = new Double_t[fNelem]; fHodoNegInvAdcAdc = new Double_t[fNelem]; fHodoSigma = new Double_t[fNelem]; + + //New Time-Walk Calibration Parameters + fHodoVelFit=new Double_t [fNelem]; + fHodoCableFit=new Double_t [fNelem]; + fHodo_LCoeff=new Double_t [fNelem]; + fHodoPos_c1=new Double_t [fNelem]; + fHodoNeg_c1=new Double_t [fNelem]; + fHodoPos_c2=new Double_t [fNelem]; + fHodoNeg_c2=new Double_t [fNelem]; + for(Int_t j=0;j<(Int_t) fNelem;j++) { Int_t index=((THcHodoscope *)GetParent())->GetScinIndex(fPlaneNum-1,j); fHodoPosAdcTimeWindowMin[j] = ((THcHodoscope*) GetParent())->GetHodoPosAdcTimeWindowMin(index); @@ -289,11 +307,25 @@ Int_t THcScintillatorPlane::ReadDatabase( const TDatime& date ) fHodoPosInvAdcAdc[j] = ((THcHodoscope *)GetParent())->GetHodoPosInvAdcAdc(index); fHodoNegInvAdcAdc[j] = ((THcHodoscope *)GetParent())->GetHodoNegInvAdcAdc(index); fHodoVelLight[j] = ((THcHodoscope *)GetParent())->GetHodoVelLight(index); + //Get Time-Walk correction param + fHodoVelFit[j] = ((THcHodoscope *)GetParent())->GetHodoVelFit(index); + fHodoCableFit[j] = ((THcHodoscope *)GetParent())->GetHodoCableFit(index); + fHodo_LCoeff[j] = ((THcHodoscope *)GetParent())->GetHodoLCoeff(index); + fHodoPos_c1[j] = ((THcHodoscope *)GetParent())->GetHodoPos_c1(index); + fHodoNeg_c1[j] = ((THcHodoscope *)GetParent())->GetHodoNeg_c1(index); + fHodoPos_c2[j] = ((THcHodoscope *)GetParent())->GetHodoPos_c2(index); + fHodoNeg_c2[j] = ((THcHodoscope *)GetParent())->GetHodoNeg_c2(index); + Double_t possigma = ((THcHodoscope *)GetParent())->GetHodoPosSigma(index); Double_t negsigma = ((THcHodoscope *)GetParent())->GetHodoNegSigma(index); fHodoSigma[j] = TMath::Sqrt(possigma*possigma+negsigma*negsigma)/2.0; + + + + } + fTdc_Thrs = ((THcHodoscope *)GetParent())->GetTDCThrs(); // cout <<" plane num = "<<fPlaneNum<<endl; // cout <<" nelem = "<<fNelem<<endl; // cout <<" zpos = "<<fZpos<<endl; @@ -343,7 +375,9 @@ Int_t THcScintillatorPlane::ReadDatabase( const TDatime& date ) fGoodNegTdcTimeCorr = vector<Double_t> (fNelem, 0.0); fGoodPosTdcTimeTOFCorr = vector<Double_t> (fNelem, 0.0); fGoodNegTdcTimeTOFCorr = vector<Double_t> (fNelem, 0.0); - + fGoodPosTdcTimeWalkCorr = vector<Double_t> (fNelem, 0.0); + fGoodNegTdcTimeWalkCorr = vector<Double_t> (fNelem, 0.0); + fGoodDiffDistTrack = vector<Double_t> (fNelem, 0.0); return kOK; } //_____________________________________________________________________________ @@ -438,10 +472,12 @@ Int_t THcScintillatorPlane::DefineVariables( EMode mode ) {"GoodNegTdcTimeUnCorr", "List of negative TDC values (passed TDC && ADC Min and Max cuts for either end)", "fGoodNegTdcTimeUnCorr"}, //Units ns {"GoodNegTdcTimeCorr", "List of negative corrected TDC values (corrected for PMT offset and ADC)", "fGoodNegTdcTimeCorr"}, {"GoodNegTdcTimeTOFCorr", "List of negative corrected TDC values (corrected for TOF)", "fGoodNegTdcTimeTOFCorr"}, + {"GoodNegTdcTimeWalkCorr", "List of negative corrected TDC values (corrected for Time-Walk)", "fGoodNegTdcTimeWalkCorr"}, {"GoodNegAdcPulseInt", "List of negative ADC values (passed TDC && ADC Min and Max cuts for either end)", "fGoodNegAdcPulseInt"}, {"GoodPosTdcTimeUnCorr", "List of positive TDC values (passed TDC && ADC Min and Max cuts for either end)", "fGoodPosTdcTimeUnCorr"}, {"GoodPosTdcTimeCorr", "List of positive corrected TDC values (corrected for PMT offset and ADC)", "fGoodPosTdcTimeCorr"}, {"GoodPosTdcTimeTOFCorr", "List of positive corrected TDC values (corrected for TOF)", "fGoodPosTdcTimeTOFCorr"}, + {"GoodPosTdcTimeWalkCorr", "List of positive corrected TDC values (corrected for Time-Walk)", "fGoodPosTdcTimeWalkCorr"}, {"GoodPosAdcPulseInt", "List of positive ADC values (passed TDC && ADC Min and Max cuts for either end)", "fGoodPosAdcPulseInt"}, {"GoodPosAdcPulseAmp", "List of positive ADC peak amp (passed TDC && ADC Min and Max cuts for either end)", "fGoodPosAdcPulseAmp"}, {"GoodNegAdcPulseAmp", "List of Negative ADC peak amp (passed TDC && ADC Min and Max cuts for either end)", "fGoodNegAdcPulseAmp"}, @@ -450,6 +486,7 @@ Int_t THcScintillatorPlane::DefineVariables( EMode mode ) {"GoodPosAdcTdcDiffTime", "List of positive TDC - ADC time (passed TDC && ADC Min and Max cuts for either end)", "fGoodPosAdcTdcDiffTime"}, {"GoodNegAdcTdcDiffTime", "List of Negative TDC - ADC time (passed TDC && ADC Min and Max cuts for either end)", "fGoodNegAdcTdcDiffTime"}, {"DiffDisTrack", "Difference between track and scintillator position (cm)", "fHitDistance"}, + {"DiffDisTrackCorr", "TW Corrected Dist Difference between track and scintillator position (cm)", "fGoodDiffDistTrack"}, {"TrackXPos", "Track X position at plane (cm)", "fTrackXPosition"}, {"TrackYPos", "Track Y position at plane (cm)", "fTrackYPosition"}, //{"ngoodhits", "Number of paddle hits (passed tof tolerance and used to determine the focal plane time )", "GetNGoodHits() "}, @@ -537,14 +574,19 @@ void THcScintillatorPlane::Clear( Option_t* ) fGoodPosTdcTimeUnCorr.at(ielem) = kBig; fGoodPosTdcTimeCorr.at(ielem) = kBig; fGoodPosTdcTimeTOFCorr.at(ielem) = kBig; + fGoodPosTdcTimeWalkCorr.at(ielem) = kBig; } for (UInt_t ielem = 0; ielem < fGoodNegTdcTimeUnCorr.size(); ielem++) { fGoodNegTdcTimeUnCorr.at(ielem) = kBig; fGoodNegTdcTimeCorr.at(ielem) = kBig; fGoodNegTdcTimeTOFCorr.at(ielem) = kBig; + fGoodNegTdcTimeWalkCorr.at(ielem) = kBig; } + for (UInt_t ielem = 0; ielem < fGoodDiffDistTrack.size(); ielem++) { + fGoodDiffDistTrack.at(ielem) = kBig; + } fpTime = -1.e4; fHitDistance = kBig; @@ -926,13 +968,27 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) ((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPosADCtime(adctime_pos); ((THcHodoHit*) fHodoHits->At(fNScinHits))->SetNegADCtime(adctime_neg); + //CALCULATE TIME-WALK CORRECTIONS HERE!!!! //Define GoodTdcUnCorrTime if(btdcraw_pos) { fGoodPosTdcTimeUnCorr.at(padnum-1) = tdc_pos*fScinTdcToTime; + //tw_corr_pos = fHodoPos_c1[padnum-1]/pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[padnum-1]) - fHodoPos_c1[padnum-1]/pow(200./fTdc_Thrs, fHodoPos_c2[padnum-1]); + + tw_corr_pos = 1./pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[padnum-1]) - 1./pow(200./fTdc_Thrs, fHodoPos_c2[padnum-1]); + + fGoodPosTdcTimeWalkCorr.at(padnum-1) = tdc_pos*fScinTdcToTime -tw_corr_pos; } if(btdcraw_neg) { fGoodNegTdcTimeUnCorr.at(padnum-1) = tdc_neg*fScinTdcToTime; + + //tw_corr_neg = fHodoNeg_c1[padnum-1]/pow(adcamp_neg/fTdc_Thrs,fHodoNeg_c2[padnum-1]) - fHodoNeg_c1[padnum-1]/pow(200./fTdc_Thrs, fHodoNeg_c2[padnum-1]); + + tw_corr_neg = 1./pow(adcamp_neg/fTdc_Thrs,fHodoNeg_c2[padnum-1]) - 1./pow(200./fTdc_Thrs, fHodoNeg_c2[padnum-1]); + + fGoodNegTdcTimeWalkCorr.at(padnum-1) = tdc_neg*fScinTdcToTime -tw_corr_neg; + } + // Do corrections if valid TDC on both ends of bar if(btdcraw_pos && btdcraw_neg) { // Do the pulse height correction to the time. (Position dependent corrections later) @@ -944,17 +1000,42 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) timec_neg = tdc_neg*fScinTdcToTime - fHodoNegInvAdcOffset[index] - fHodoNegInvAdcAdc[index]/TMath::Sqrt(TMath::Max(20.0*.020,adcint_neg)); - } else { // Old style - timec_pos = tdc_pos*fScinTdcToTime - fHodoPosPhcCoeff[index]* - TMath::Sqrt(TMath::Max(0.0,adcint_pos/fHodoPosMinPh[index]-1.0)) - - fHodoPosTimeOffset[index]; - timec_neg = tdc_neg*fScinTdcToTime - fHodoNegPhcCoeff[index]* - TMath::Sqrt(TMath::Max(0.0,adcint_neg/fHodoNegMinPh[index]-1.0)) - - fHodoNegTimeOffset[index]; + } else { // FADC style + timec_pos = tdc_pos*fScinTdcToTime -tw_corr_pos + fHodo_LCoeff[index]; + timec_neg = tdc_neg*fScinTdcToTime -tw_corr_neg- 2*fHodoCableFit[index] + fHodo_LCoeff[index]; } + + Double_t TWCorrDiff = fGoodNegTdcTimeWalkCorr.at(padnum-1) - 2*fHodoCableFit[index] - fGoodPosTdcTimeWalkCorr.at(padnum-1); + + Double_t fHitDistCorr = 0.5*TWCorrDiff*fHodoVelFit[index]; + + /*Debug*/ + //cout << "*****************" << endl; + //cout << "fPlNum: " << fPlaneNum << endl; + //cout << "*****************" << endl; + //cout << "Paddle index: " << index << endl; + //cout << "pos_sigma: " << fHodoSigma[index] << endl; + //cout << "fZPos: " << fZpos << endl; + //cout << "fDzPos: " << fDzpos << endl; + //cout << "Zcorr = fZpos+(index%2)*fDzpos = " << fZpos+(index%2)*fDzpos << endl; + + //cout << Form("****fHodo_LCoeff[%d]", index) << fHodo_LCoeff[index] << endl; + //cout << Form("****fHodoCableFit[%d]", index) << fHodoCableFit[index] << endl; + //cout << Form("****fHodoVelFit[%d]", index) << fHodoVelFit[index] << endl; + //cout << Form("****c1_Pos/Neg[%d]", index) << fHodoPos_c1[index] << " / " << fHodoNeg_c1[index] << endl; + //cout << Form("****c2_Pos/Neg[%d]", index) << fHodoPos_c2[index] << " / " << fHodoNeg_c2[index] << endl; + //cout << "TW Corr val. +/-: " << tw_corr_pos << " / " << tw_corr_neg << endl; + //cout << "TW UnCorr +/-: " << fGoodPosTdcTimeUnCorr.at(padnum-1) << " / " << fGoodNegTdcTimeUnCorr.at(padnum-1) << endl; + //cout << "TW Corr +/-: " << fGoodPosTdcTimeWalkCorr.at(padnum-1) << " / " << fGoodNegTdcTimeWalkCorr.at(padnum-1) << endl; + + + fGoodDiffDistTrack.at(index) = fHitDistCorr; // Find hit position using ADCs // If postime larger, then hit was nearer negative side. - Double_t vellight=fHodoVelLight[index]; + + Double_t vellight=fHodoVelLight[index]; //read from hodo_cuts.param, where it is set fixed to 15.0 + //Double_t vellight=fHodoVelFit[index]; //use scin prop vel. values from hodo_calibVp_run#.param file + Double_t dist_from_center=0.5*(timec_neg-timec_pos)*vellight; Double_t scint_center=0.5*(fPosLeft+fPosRight); Double_t hit_position=scint_center+dist_from_center; @@ -975,9 +1056,9 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) negtime = timec_neg - (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal); } } else { - postime=timec_pos-(fPosLeft-hit_position)/fHodoVelLight[index]; - negtime=timec_neg-(hit_position-fPosRight)/fHodoVelLight[index]; - scin_corrected_time = 0.5*(postime+negtime); + scin_corrected_time = 0.5*(timec_neg+timec_pos); // add constants for each paddle, 25ns, 25 + zpos, . . . //remove propagation time + timec_pos= scin_corrected_time; + timec_neg= scin_corrected_time; if (fCosmicFlag) { postime = timec_pos + (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal); negtime = timec_neg + (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal); @@ -990,6 +1071,8 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) ((THcHodoHit*) fHodoHits->At(fNScinHits))->SetCorrectedTimes(timec_pos,timec_neg, postime, negtime, scin_corrected_time); + ((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPosADCpeak(adcamp_pos); // need for new TWCOrr + ((THcHodoHit*) fHodoHits->At(fNScinHits))->SetNegADCpeak(adcamp_neg); // need for new TWCOrr fGoodPosTdcTimeCorr.at(padnum-1) = timec_pos; fGoodNegTdcTimeCorr.at(padnum-1) = timec_neg; fGoodPosTdcTimeTOFCorr.at(padnum-1) = postime; @@ -1003,10 +1086,8 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) timec_pos = tdc_pos*fScinTdcToTime - fHodoPosInvAdcOffset[index] - fHodoPosInvAdcAdc[index]/TMath::Sqrt(TMath::Max(20.0*.020,adcint_pos)); - } else { // Old style - timec_pos = tdc_pos*fScinTdcToTime - fHodoPosPhcCoeff[index]* - TMath::Sqrt(TMath::Max(0.0,adcint_pos/fHodoPosMinPh[index]-1.0)) - - fHodoPosTimeOffset[index]; + } else { // FADC style + timec_pos = tdc_pos*fScinTdcToTime -tw_corr_pos + fHodo_LCoeff[index]; } } if(btdcraw_neg) { @@ -1014,15 +1095,15 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) timec_neg = tdc_neg*fScinTdcToTime - fHodoNegInvAdcOffset[index] - fHodoNegInvAdcAdc[index]/TMath::Sqrt(TMath::Max(20.0*.020,adcint_neg)); - } else { // Old style - timec_neg = tdc_neg*fScinTdcToTime - fHodoNegPhcCoeff[index]* - TMath::Sqrt(TMath::Max(0.0,adcint_neg/fHodoNegMinPh[index]-1.0)) - - fHodoNegTimeOffset[index]; + } else { // FADC style + timec_neg = tdc_neg*fScinTdcToTime -tw_corr_neg- 2*fHodoCableFit[index] + fHodo_LCoeff[index]; } } ((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPaddleCenter(fPosCenter[index]); ((THcHodoHit*) fHodoHits->At(fNScinHits))->SetCorrectedTimes(timec_pos,timec_neg, timec_pos,timec_neg,0.0); + ((THcHodoHit*) fHodoHits->At(fNScinHits))->SetNegADCpeak(adcamp_neg); // needed for new TWCOrr + ((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPosADCpeak(adcamp_pos); // needed for new TWCOrr fGoodPosTdcTimeCorr.at(padnum-1) = timec_pos; fGoodNegTdcTimeCorr.at(padnum-1) = timec_neg; fGoodPosTdcTimeTOFCorr.at(padnum-1) = timec_pos; diff --git a/src/THcScintillatorPlane.h b/src/THcScintillatorPlane.h index e7b38fde4eb6973f55e40e565b1da5d8b5dfe521..e39d514712ab4c826542a860e03920fd23976006 100644 --- a/src/THcScintillatorPlane.h +++ b/src/THcScintillatorPlane.h @@ -162,7 +162,10 @@ class THcScintillatorPlane : public THaSubDetector { vector<Double_t> fGoodNegTdcTimeCorr; vector<Double_t> fGoodNegTdcTimeTOFCorr; - + //Time Walk Corrected + vector<Double_t> fGoodPosTdcTimeWalkCorr; + vector<Double_t> fGoodNegTdcTimeWalkCorr; + vector<Double_t> fGoodDiffDistTrack; Int_t fDebugAdc; Double_t fHitDistance; @@ -222,6 +225,19 @@ class THcScintillatorPlane : public THaSubDetector { Double_t *fHodoNegInvAdcLinear; Double_t *fHodoPosInvAdcAdc; Double_t *fHodoNegInvAdcAdc; + //Time-Walk Parameters + Double_t* fHodoVelFit; + Double_t* fHodoCableFit; + Double_t* fHodo_LCoeff; + Double_t* fHodoPos_c1; + Double_t* fHodoNeg_c1; + Double_t* fHodoPos_c2; + Double_t* fHodoNeg_c2; + Double_t fTdc_Thrs; + + Double_t tw_corr_pos; + Double_t tw_corr_neg; + Double_t *fHodoSigma; Double_t fTolerance; /* need this for Focal Plane Time estimation */ Double_t fFptime;