diff --git a/podd b/podd index 993f6aaf52920ee7c8f441c88be19343823bcb38..afce4df317bf88ff161519ee1ee94bd9eb85744f 160000 --- a/podd +++ b/podd @@ -1 +1 @@ -Subproject commit 993f6aaf52920ee7c8f441c88be19343823bcb38 +Subproject commit afce4df317bf88ff161519ee1ee94bd9eb85744f diff --git a/src/THcCherenkov.cxx b/src/THcCherenkov.cxx index 31a3418027ffead7028a3198e1345e62a380520d..da3232e95a41031f7f06dd6a11f5200940aa3162 100644 --- a/src/THcCherenkov.cxx +++ b/src/THcCherenkov.cxx @@ -168,7 +168,8 @@ Int_t THcCherenkov::ReadDatabase( const TDatime& date ) fNelem = (Int_t)gHcParms->Find(parname)->GetValue(); // class. // fNelem = 2; // Default if not defined - + fCerNRegions = 3; + fNPMT = new Int_t[fNelem]; fADC = new Double_t[fNelem]; fADC_P = new Double_t[fNelem]; @@ -188,7 +189,6 @@ Int_t THcCherenkov::ReadDatabase( const TDatime& date ) fPedLimit = new Int_t[fNelem]; fPedMean = new Double_t[fNelem]; - fCerNRegions = 3; // This value should be in parameter file fCerTrackCounter = new Int_t [fCerNRegions]; fCerFiredCounter = new Int_t [fCerNRegions]; @@ -197,22 +197,22 @@ Int_t THcCherenkov::ReadDatabase( const TDatime& date ) fCerFiredCounter[ireg] = 0; } - fCerRegionsValueMax = fCerNRegions * 8; // This value 8 should also be in paramter file fCerRegionValue = new Double_t [fCerRegionsValueMax]; DBRequest list[]={ - {"cer_adc_to_npe", fGain, kDouble, (UInt_t) fNelem}, // Ahmed - {"cer_ped_limit", fPedLimit, kInt, (UInt_t) fNelem}, // Ahmed - {"cer_width", fCerWidth, kDouble, (UInt_t) fNelem}, // Ahmed - {"cer_chi2max", &fCerChi2Max, kDouble}, // Ahmed - {"cer_beta_min", &fCerBetaMin, kDouble}, // Ahmed - {"cer_beta_max", &fCerBetaMax, kDouble}, // Ahmed - {"cer_et_min", &fCerETMin, kDouble}, // Ahmed - {"cer_et_max", &fCerETMax, kDouble}, // Ahmed - {"cer_mirror_zpos", &fCerMirrorZPos, kDouble}, // Ahmed - {"cer_region", &fCerRegionValue[0], kDouble, (UInt_t) fCerRegionsValueMax}, // Ahmed - {"cer_threshold", &fCerThresh, kDouble}, // Ahmed + {"cer_adc_to_npe", fGain, kDouble, (UInt_t) fNelem}, + {"cer_ped_limit", fPedLimit, kInt, (UInt_t) fNelem}, + {"cer_width", fCerWidth, kDouble, (UInt_t) fNelem}, + {"cer_chi2max", &fCerChi2Max, kDouble}, + {"cer_beta_min", &fCerBetaMin, kDouble}, + {"cer_beta_max", &fCerBetaMax, kDouble}, + {"cer_et_min", &fCerETMin, kDouble}, + {"cer_et_max", &fCerETMax, kDouble}, + {"cer_mirror_zpos", &fCerMirrorZPos, kDouble}, + {"cer_region", &fCerRegionValue[0], kDouble, (UInt_t) fCerRegionsValueMax}, + {"cer_threshold", &fCerThresh, kDouble}, + // {"cer_regions", &fCerNRegions, kInt}, {0} }; @@ -220,6 +220,7 @@ Int_t THcCherenkov::ReadDatabase( const TDatime& date ) fIsInit = true; + for (Int_t i1 = 0; i1 < fCerNRegions; i1++ ) { cout << "Region " << i1 << endl; for (Int_t i2 = 0; i2 < 8; i2++ ) { @@ -250,12 +251,12 @@ Int_t THcCherenkov::DefineVariables( EMode mode ) // No. They show up in tree as Ndata.H.aero.postdchits for example RVarDef vars[] = { - {"phototubes", "Nuber of Cherenkov photo tubes", "fNPMT"}, - {"adc", "Raw ADC values", "fADC"}, - {"adc_p", "Pedestal Subtracted ADC values", "fADC_P"}, - {"npe", "Number of Photo electrons", "fNPE"}, - {"npesum", "Sum of Number of Photo electrons", "fNPEsum"}, - {"ncherhit", "Number of Hits(Cherenkov)", "fNCherHit"}, + {"phototubes", "Nuber of Cherenkov photo tubes", "fNPMT"}, + {"adc", "Raw ADC values", "fADC"}, + {"adc_p", "Pedestal Subtracted ADC values", "fADC_P"}, + {"npe", "Number of Photo electrons", "fNPE"}, + {"npesum", "Sum of Number of Photo electrons", "fNPEsum"}, + {"ncherhit", "Number of Hits(Cherenkov)", "fNCherHit"}, {"certrackcounter", "Tracks inside Cherenkov region", "fCerTrackCounter"}, {"cerfiredcounter", "Tracks with engough Cherenkov NPEs ", "fCerFiredCounter"}, { 0 } @@ -374,7 +375,7 @@ Int_t THcCherenkov::FineProcess( TClonesArray& tracks ) THaTrack* theTrack = dynamic_cast<THaTrack*>( tracks.At(0) ); if (!theTrack) return -1; - + if ( ( ( tracks.GetLast() + 1 ) == 1 ) && ( theTrack->GetChi2()/theTrack->GetNDoF() > 0. ) && ( theTrack->GetChi2()/theTrack->GetNDoF() < fCerChi2Max ) && @@ -405,18 +406,12 @@ Int_t THcCherenkov::FineProcess( TClonesArray& tracks ) fCerTrackCounter[ir] ++; // * increment the 'did fire' counters - if ( fNPEsum > fCerThresh ) { fCerFiredCounter[ir] ++; } - } - - } // loop over regions - // cout << endl; - + } // loop over regions } - } return 0; diff --git a/src/THcDC.cxx b/src/THcDC.cxx index f1e8290bdd9c7ea80e97b3b4fe731ece363a055a..d714583e0b63a73e83acdb681426e45a19b4a42e 100644 --- a/src/THcDC.cxx +++ b/src/THcDC.cxx @@ -346,6 +346,7 @@ Int_t THcDC::DefineVariables( EMode mode ) // Register variables in global list RVarDef vars[] = { + { "stubtest", "stub test", "fStubTest" }, { "nhit", "Number of DC hits", "fNhits" }, { "tnhit", "Number of good DC hits", "fNthits" }, { "trawhit", "Number of true raw DC hits", "fN_True_RawHits" }, @@ -422,6 +423,7 @@ inline void THcDC::ClearEvent() { // Reset per-event data. + fStubTest = 0; fNhits = 0; fNthits = 0; fN_True_RawHits=0; @@ -695,6 +697,7 @@ void THcDC::LinkStubs() && (TMath::Abs(dposyp) < fYptTrCriterion)) { if(newtrack) { assert(sptracks==0); + fStubTest = 1; //stubtest=1; Used in h_track_tests.f // Make a new track if there are not to many if(fNDCTracks < MAXTRACKS) { diff --git a/src/THcDC.h b/src/THcDC.h index 71bd4963ff9de5849de4081579ed0d7bb79ff54b..d11f2333ce479486419f285fc014a2d062889c96 100644 --- a/src/THcDC.h +++ b/src/THcDC.h @@ -103,6 +103,7 @@ protected: // Was used for SOS in ENGINE. // Per-event data + Int_t fStubTest; Int_t fNhits; Int_t fNthits; Int_t fN_True_RawHits; diff --git a/src/THcDriftChamber.cxx b/src/THcDriftChamber.cxx index b9ed691797b48154cf905f7ca2d4330b3d7bed22..95906b6184edb64589a5dac2464e2abf9c787430 100644 --- a/src/THcDriftChamber.cxx +++ b/src/THcDriftChamber.cxx @@ -224,6 +224,8 @@ Int_t THcDriftChamber::DefineVariables( EMode mode ) // Register variables in global list RVarDef vars[] = { + { "maxhits", "Maximum hits allowed", "fMaxHits" }, + { "spacepoints", "Space points of DC", "fNSpacePoints" }, { "nhit", "Number of DC hits", "fNhits" }, { "trawhit", "Number of True Raw hits", "fN_True_RawHits" }, { 0 } diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx index 16c30403e4f8a637a2deea47cee76d39fecf6be8..6f821785efca275c73d80e861859e0f1e3169f19 100644 --- a/src/THcHodoscope.cxx +++ b/src/THcHodoscope.cxx @@ -259,8 +259,11 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date ) // Int_t plen=strlen(parname); cout << " readdatabse hodo fnplanes = " << fNPlanes << endl; + fBetaP = 0.; + fBetaNoTrk = 0.; + fBetaNoTrkChiSq = 0.; + fNPaddle = new UInt_t [fNPlanes]; - fFPTime = new Double_t [fNPlanes]; fPlaneCenter = new Double_t[fNPlanes]; fPlaneSpacing = new Double_t[fNPlanes]; @@ -430,13 +433,16 @@ Int_t THcHodoscope::DefineVariables( EMode mode ) RVarDef vars[] = { // Move these into THcHallCSpectrometer using track fTracks - {"fpHitsTime", "Time at focal plane from all hits", "fFPTime"}, - {"starttime", "Hodoscope Start Time", "fStartTime"}, - {"goodstarttime", "Hodoscope Good Start Time", "fGoodStartTime"}, - {"goodscinhit", "Hit in fid area", "fGoodScinHits"}, - // {"goodscinhitx", "Hit in fid x range", "fGoodScinHitsX"}, - {"scinshould", "Total scin Hits in fid area", "fScinShould"}, - {"scindid", "Total scin Hits in fid area with a track", "fScinDid"}, + {"betap", "betaP", "fBetaP"}, + {"betanotrack", "Beta from scintillator hits", "fBetaNoTrk"}, + {"betachisqnotrack", "Chi square of beta from scintillator hits", "fBetaNoTrkChiSq"}, + {"fpHitsTime", "Time at focal plane from all hits", "fFPTime"}, + {"starttime", "Hodoscope Start Time", "fStartTime"}, + {"goodstarttime", "Hodoscope Good Start Time", "fGoodStartTime"}, + {"goodscinhit", "Hit in fid area", "fGoodScinHits"}, + // {"goodscinhitx", "Hit in fid x range", "fGoodScinHitsX"}, + {"scinshould", "Total scin Hits in fid area", "fScinShould"}, + {"scindid", "Total scin Hits in fid area with a track", "fScinDid"}, { 0 } }; return DefineVarsFromList( vars, mode ); @@ -504,12 +510,10 @@ void THcHodoscope::DeleteArrays() inline void THcHodoscope::ClearEvent() { - // Reset per-event data. - // for ( Int_t imaxhit = 0; imaxhit < MAXHODHITS; imaxhit++ ){ - // fBeta[imaxhit] = 0.; - // fBetaChisq[imaxhit] = 0.; - // } + fBetaP = 0.; + fBetaNoTrk = 0.0; + fBetaNoTrkChiSq = 0.0; for(Int_t ip=0;ip<fNPlanes;ip++) { fPlanes[ip]->Clear(); @@ -590,6 +594,7 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata ) //_____________________________________________________________________________ void THcHodoscope::EstimateFocalPlaneTime( void ) { + Int_t timehist[200]; for (Int_t i=0;i<200;i++) { @@ -629,17 +634,27 @@ void THcHodoscope::EstimateFocalPlaneTime( void ) ihit = 0; Double_t fpTimeSum = 0.0; + Int_t jhit = 0; fNfptimes=0; + fNoTrkPlaneInfo.clear(); + fNoTrkHitInfo.clear(); for(Int_t ip=0;ip<fNPlanes;ip++) { + fNoTrkPlaneInfo.push_back(NoTrkPlaneInfo()); + fNoTrkPlaneInfo[ip].goodplanetime = kFALSE; Int_t nphits=fPlanes[ip]->GetNScinHits(); TClonesArray* hodoHits = fPlanes[ip]->GetHits(); for(Int_t i=0;i<nphits;i++) { + fNoTrkHitInfo.push_back(NoTrkHitInfo()); + fNoTrkHitInfo[jhit].goodtwotimes = kFALSE; + fNoTrkHitInfo[jhit].goodscintime = kFALSE; Double_t tmin = 0.5*binmax; Double_t postime=((THcHodoHit*) hodoHits->At(i))->GetPosTOFCorrectedTime(); Double_t negtime=((THcHodoHit*) hodoHits->At(i))->GetNegTOFCorrectedTime(); if ((postime>tmin) && (postime<tmin+fTofTolerance) && (negtime>tmin) && (negtime<tmin+fTofTolerance)) { + fNoTrkHitInfo[jhit].goodtwotimes = kTRUE; + fNoTrkHitInfo[jhit].goodscintime = kTRUE; // Both tubes fired Int_t index=((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber()-1; Double_t fptime = ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime() @@ -649,11 +664,14 @@ void THcHodoscope::EstimateFocalPlaneTime( void ) // Should also fill the all FP times histogram fpTimeSum += fptime; fNfptimes++; + fNoTrkPlaneInfo[ip].goodplanetime = kTRUE; } } + jhit++; } ihit++; } + if(fNfptimes>0) { fStartTime = fpTimeSum/fNfptimes; fGoodStartTime=kTRUE; @@ -661,6 +679,89 @@ void THcHodoscope::EstimateFocalPlaneTime( void ) fStartTime = fStartTimeCenter; fGoodStartTime=kFALSE; } + + + if ( ( fNoTrkPlaneInfo[0].goodplanetime || fNoTrkPlaneInfo[1].goodplanetime ) && + ( fNoTrkPlaneInfo[2].goodplanetime || fNoTrkPlaneInfo[3].goodplanetime ) ){ + + Double_t sumW = 0.; + Double_t sumT = 0.; + Double_t sumZ = 0.; + Double_t sumZZ = 0.; + Double_t sumTZ = 0.; + Int_t ihhit = 0; + + for(Int_t ip=0;ip<fNPlanes;ip++) { + Int_t nphits=fPlanes[ip]->GetNScinHits(); + TClonesArray* hodoHits = fPlanes[ip]->GetHits(); + + for(Int_t i=0;i<nphits;i++) { + Int_t index=((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber()-1; + + if ( fNoTrkHitInfo[ihhit].goodscintime ) { + + Double_t sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) + + TMath::Power( fHodoNegSigma[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; + + sumW += scinWeight; + sumT += scinWeight * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime(); + sumZ += scinWeight * zPosition; + sumZZ += scinWeight * ( zPosition * zPosition ); + sumTZ += scinWeight * zPosition * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime(); + + } // condition of good scin time + ihhit ++; + } // loop over hits of plane + } // loop over planes + + Double_t tmp = sumW * sumZZ - sumZ * sumZ ; + Double_t t0 = ( sumT * sumZZ - sumZ * sumTZ ) / tmp ; + Double_t tmpDenom = sumW * sumTZ - sumZ * sumT; + + if ( TMath::Abs( tmpDenom ) > ( 1 / 10000000000.0 ) ) { + + fBetaNoTrk = tmp / tmpDenom; + fBetaNoTrkChiSq = 0.; + ihhit = 0; + + for (Int_t ip = 0; ip < fNPlanes; ip++ ){ // Loop over planes + Int_t nphits=fPlanes[ip]->GetNScinHits(); + TClonesArray* hodoHits = fPlanes[ip]->GetHits(); + + for(Int_t i=0;i<nphits;i++) { + Int_t index=((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber()-1; + + if ( fNoTrkHitInfo[ihhit].goodscintime ) { + + 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) ) ); + fBetaNoTrkChiSq += ( ( zPosition / fBetaNoTrk - timeDif ) * + ( zPosition / fBetaNoTrk - timeDif ) ) / ( sigma * sigma ); + + + } // condition for good scin time + ihhit++; + } // loop over hits of a plane + } // loop over planes + + Double_t pathNorm = 1.0; + + fBetaNoTrk = fBetaNoTrk * pathNorm; + fBetaNoTrk = fBetaNoTrk / 29.979; // velocity / c + + } // condition for fTmpDenom + else { + fBetaNoTrk = 0.; + fBetaNoTrkChiSq = -2.; + } // else condition for fTmpDenom + } + } //_____________________________________________________________________________ Int_t THcHodoscope::ApplyCorrections( void ) @@ -722,7 +823,7 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks ) Double_t sumFPTime = 0.; // Line 138 fNScinHit.push_back(0); Double_t p = theTrack->GetP(); // Line 142 - Double_t betaP = p/( TMath::Sqrt( p * p + fPartMass * fPartMass) ); + fBetaP = p/( TMath::Sqrt( p * p + fPartMass * fPartMass) ); //! Calculate all corrected hit times and histogram //! This uses a copy of code below. Results are save in time_pos,neg @@ -797,7 +898,7 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks ) Double_t pathp = fPlanes[ip]->GetPosLeft() - scinLongCoord; Double_t timep = ((THcHodoHit*)hodoHits->At(iphit))->GetPosCorrectedTime(); timep = timep - ( pathp / fHodoVelLight[fPIndex] ) - ( fPlanes[ip]->GetZpos() + - ( paddle % 2 ) * fPlanes[ip]->GetDzpos() ) / ( 29.979 * betaP ) * + ( paddle % 2 ) * fPlanes[ip]->GetDzpos() ) / ( 29.979 * fBetaP ) * TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() + theTrack->GetPhi() * theTrack->GetPhi() ); fTOFPInfo[iphit].time_pos = timep; @@ -811,7 +912,7 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks ) Double_t pathn = scinLongCoord - fPlanes[ip]->GetPosRight(); Double_t timen = ((THcHodoHit*)hodoHits->At(iphit))->GetNegCorrectedTime(); timen = timen - ( pathn / fHodoVelLight[fPIndex] ) - ( fPlanes[ip]->GetZpos() + - ( paddle % 2 ) * fPlanes[ip]->GetDzpos() ) / ( 29.979 * betaP ) * + ( paddle % 2 ) * fPlanes[ip]->GetDzpos() ) / ( 29.979 * fBetaP ) * TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() + theTrack->GetPhi() * theTrack->GetPhi() ); fTOFPInfo[iphit].time_neg = timen; @@ -954,7 +1055,7 @@ Int_t THcHodoscope::FineProcess( TClonesArray& tracks ) // scin_time_fp doesn't need to be an array Double_t scin_time_fp = fTOFCalc[ihhit].scin_time - ( fPlanes[ip]->GetZpos() + ( paddle % 2 ) * fPlanes[ip]->GetDzpos() ) / - ( 29.979 * betaP ) * + ( 29.979 * fBetaP ) * TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() + theTrack->GetPhi() * theTrack->GetPhi() ); diff --git a/src/THcHodoscope.h b/src/THcHodoscope.h index efc502bc3f34dd040ce714972d4ebef5f8ddebf9..d113af1fcaf15314d3176930262a0d52b78851c7 100644 --- a/src/THcHodoscope.h +++ b/src/THcHodoscope.h @@ -69,7 +69,7 @@ public: Double_t GetHodoVelLight(Int_t iii) const {return fHodoVelLight[iii];} Double_t GetStartTimeCenter() const {return fStartTimeCenter;} Double_t GetStartTimeSlop() const {return fStartTimeSlop;} - Double_t GetBetaNotrk() const {return fBetaNotrk;} + Double_t GetBetaNotrk() const {return fBetaNoTrk;} Int_t GetGoodRawPad(Int_t iii){return fTOFCalc[iii].good_raw_pad;} Int_t GetGoodRawPlane(Int_t iii){return fTOFCalc[iii].pindex;} @@ -107,7 +107,10 @@ protected: Double_t fStartTime; Int_t fNfptimes; - Double_t fBetaNotrk; + Double_t fBetaP; + + Double_t fBetaNoTrk; + Double_t fBetaNoTrkChiSq; // Per-event data // Potential Hall C parameters. Mostly here for demonstration @@ -223,6 +226,22 @@ protected: // Double_t* gain; // } fDataDest[NDEST]; // Lookup table for decoder + // Inforamtion for each plane + struct NoTrkPlaneInfo { + Bool_t goodplanetime; + NoTrkPlaneInfo () : goodplanetime(kFALSE) {} + }; + std::vector<NoTrkPlaneInfo> fNoTrkPlaneInfo; + + // Inforamtion for each plane + struct NoTrkHitInfo { + Bool_t goodtwotimes; + Bool_t goodscintime; + NoTrkHitInfo () : goodtwotimes(kFALSE) {} + }; + std::vector<NoTrkHitInfo> fNoTrkHitInfo; + + // Used in TOF calculation (FineProcess) to hold information about hits // within a given plane struct TOFPInfo { diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx index f56ef154a3058ac2bf1fbc462eb28ab9d64fdefc..5a6ada39b1ac4b742cb50011945ffa71aa246abb 100644 --- a/src/THcScintillatorPlane.cxx +++ b/src/THcScintillatorPlane.cxx @@ -222,16 +222,12 @@ Int_t THcScintillatorPlane::DefineVariables( EMode mode ) // Register variables in global list RVarDef vars[] = { - {"postdchits", "List of Positive TDC hits", - "frPosTDCHits.THcSignalHit.GetPaddleNumber()"}, - {"negtdchits", "List of Negative TDC hits", - "frNegTDCHits.THcSignalHit.GetPaddleNumber()"}, - {"posadchits", "List of Positive ADC hits", - "frPosADCHits.THcSignalHit.GetPaddleNumber()"}, - {"negadchits", "List of Negative ADC hits", - "frNegADCHits.THcSignalHit.GetPaddleNumber()"}, - // {"fptime", "Time at focal plane", - // "GetFpTime()"}, + {"postdchits", "List of Positive TDC hits", "frPosTDCHits.THcSignalHit.GetPaddleNumber()"}, + {"negtdchits", "List of Negative TDC hits", "frNegTDCHits.THcSignalHit.GetPaddleNumber()"}, + {"posadchits", "List of Positive ADC hits", "frPosADCHits.THcSignalHit.GetPaddleNumber()"}, + {"negadchits", "List of Negative ADC hits", "frNegADCHits.THcSignalHit.GetPaddleNumber()"}, + // {"fptime", "Time at focal plane", // "GetFpTime()"}, + {"nhits", "Number of hits", "GetNScinHits() "}, { 0 } };