diff --git a/podd b/podd index ac15a305d5c3f7e05ce24ef47e77222327774e34..a7ebab4863214012db9ef82561c7b68ec4b1a4ea 160000 --- a/podd +++ b/podd @@ -1 +1 @@ -Subproject commit ac15a305d5c3f7e05ce24ef47e77222327774e34 +Subproject commit a7ebab4863214012db9ef82561c7b68ec4b1a4ea diff --git a/src/THcAerogel.cxx b/src/THcAerogel.cxx index 41152977a1b1dcf3ddfdac9f1c306d2f23091f39..5dba21399921887959986dcf4bcd95e5900bd831 100644 --- a/src/THcAerogel.cxx +++ b/src/THcAerogel.cxx @@ -113,8 +113,8 @@ Int_t THcAerogel::ReadDatabase( const TDatime& date ) {"aero_neg_gain", fPosGain, kDouble}, {"aero_pos_ped_limit", fPosPedLimit, kInt}, {"aero_neg_ped_limit", fNegPedLimit, kInt}, - {"aero_pos_ped_mean", fPosPedMean, kDouble}, - {"aero_neg_ped_mean", fNegPedMean, kDouble}, + // {"aero_pos_ped_mean", fPosPedMean, kDouble}, + // {"aero_neg_ped_mean", fNegPedMean, kDouble}, {0} }; gHcParms->LoadParmValues((DBRequest*)&list,prefix); @@ -155,7 +155,7 @@ Int_t THcAerogel::DefineVariables( EMode mode ) {"neg_npe","PEs PE Negative Tube","fNegNpe"}, {"pos_npe_sum", "Total Positive Tube PEs", "fPosNpeSum"}, {"neg_npe_sum", "Total Negative Tube PEs", "fNegNpeSum"}, - {"npe_sum", "Total PEs", ""}, + {"npe_sum", "Total PEs", "fNpeSum"}, {"ntdc_pos_hits", "Number of Positive Tube Hits", "fNTDCPosHits"}, {"ntdc_neg_hits", "Number of Negative Tube Hits", "fNTDCNegHits"}, {"ngood_hits", "Total number of good hits", "fNGoodHits"}, diff --git a/src/THcDriftChamberPlane.cxx b/src/THcDriftChamberPlane.cxx index f836582ce22e81fbc21117a9923e256b2b16741a..cae65633687f274352318e6f5f9ac3d5214d4a7a 100644 --- a/src/THcDriftChamberPlane.cxx +++ b/src/THcDriftChamberPlane.cxx @@ -188,7 +188,7 @@ Int_t THcDriftChamberPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) // Would be nice to have a way to determine that the hodoscope decode was // actually called for this event. if( fglHod ) StartTime = fglHod->GetStartTime(); - cout << "Start time " << StartTime << endl; + //cout << "Start time " << StartTime << endl; //Int_t nTDCHits=0; fHits->Clear(); diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx index c5c1fd68cf1dbcafd2865f64cec9f09046af7977..de6278f03603d973137bbfc7a5e52bb3ecf8e474 100644 --- a/src/THcHodoscope.cxx +++ b/src/THcHodoscope.cxx @@ -43,7 +43,8 @@ THcHodoscope::THcHodoscope( const char* name, const char* description, //fTrackProj = new TClonesArray( "THaTrackProj", 5 ); // Construct the planes fNPlanes = 0; // No planes until we make them - + fStartTime=-1e5; + fGoodStartTime=kFALSE; } //_____________________________________________________________________________ @@ -416,7 +417,6 @@ Int_t THcHodoscope::ReadDatabase( const TDatime& date ) } cout <<endl; } - cout <<endl<<endl; // check fHodoPosPhcCoeff /* @@ -592,16 +592,30 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata ) // Let each plane get its hits Int_t nexthit = 0; + Int_t nfptimes=0; + for(Int_t ip=0;ip<fNPlanes;ip++) { // nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit); // GN: select only events that have reasonable TDC values to start with // as per the Engine h_strip_scin.f nexthit = fPlanes[ip]->ProcessHits(fRawHitList,nexthit); - if (fPlanes[ip]->GetNScinHits()>0) { + if (fPlanes[ip]->GetNScinHits()>0) { fPlanes[ip]->PulseHeightCorrection(); + ///cout <<"Plane "<<ip<<" number of fpTimes = "<<fPlanes[ip]->GetFpTimeHits()<<endl; + for (Int_t ifptimes=0;ifptimes<fPlanes[ip]->GetFpTimeHits();ifptimes++) { + fStartTime=fStartTime+fPlanes[ip]->GetFpTime(ifptimes); + nfptimes++; } + } } - + if (nfptimes>0) { + fStartTime=fStartTime/nfptimes; + fGoodStartTime=kTRUE; + } else { + fGoodStartTime=kFALSE; + fStartTime=fStartTimeCenter; + } + cout <<" stats = "<<fGoodStartTime<<" "<<nfptimes<<" fStartTime = "<<fStartTime<<endl; // fRawHitList is TClones array of THcHodoscopeHit objects #if 0 for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) { @@ -613,7 +627,7 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata ) cout << endl; #endif - fStartTime = 500; // Drift Chamber will need this + /// fStartTime = 500; // Drift Chamber will need this return nhits; } diff --git a/src/THcHodoscope.h b/src/THcHodoscope.h index e17cef190cd85459f30fe3d2f4e02c9b04f3dc85..7d18bc8d808b8ea6cead145c8c615e79544b0504 100644 --- a/src/THcHodoscope.h +++ b/src/THcHodoscope.h @@ -30,8 +30,8 @@ public: virtual Int_t FineProcess( TClonesArray& tracks ); virtual Int_t ApplyCorrections( void ); - // Int_t GetNHits() const { return fNhit; } Double_t GetStartTime() const { return fStartTime; } + Bool_t IsStartTimeGood() const {return fGoodStartTime;}; Int_t GetScinIndex(Int_t nPlane, Int_t nPaddle); Int_t GetScinIndex(Int_t nSide, Int_t nPlane, Int_t nPaddle); Double_t GetPathLengthCentral(); @@ -60,7 +60,7 @@ protected: // Calibration // Per-event data - + Bool_t fGoodStartTime; Double_t fStartTime; // Per-event data @@ -70,7 +70,7 @@ protected: Double_t fStartTimeCenter, fStartTimeSlop, fScinTdcToTime; Double_t fTofTolerance; Double_t fPathLengthCentral; - Int_t fScinTdcMin, fScinTdcMax; // min and max TDC values + Double_t fScinTdcMin, fScinTdcMax; // min and max TDC values char** fPlaneNames; Int_t* fNPaddle; // Number of paddles per plane diff --git a/src/THcParmList.cxx b/src/THcParmList.cxx index ac4c119dac1171c5e2a1a178ad8e3a3e07ca41f1..0d9a845a00eaabcafaaeb5d261fb7c274d2d161a 100644 --- a/src/THcParmList.cxx +++ b/src/THcParmList.cxx @@ -16,6 +16,8 @@ #include "THaVar.h" #include "THaFormula.h" +#include "TMath.h" + /* #incluce <algorithm> include <fstream> include <cstring> */ #include <iostream> @@ -430,43 +432,53 @@ Int_t THcParmList::LoadParmValues(const DBRequest* list, const char* prefix) string keystr(prefix); keystr.append(ti->name); const char* key = keystr.c_str(); /// cout <<"Now at "<<ti->name<<endl; - if (ti->nelem>1) { - // it is an array, use the appropriateinterface - switch (ti->type) { - case (kDouble) : - this_cnt = GetArray(key,static_cast<Double_t*>(ti->var),ti->nelem); - break; - case (kInt) : - this_cnt = GetArray(key,static_cast<Int_t*>(ti->var),ti->nelem); - break; - default: - Error("THcParmList","Invalid type to read %s",key); - break; - } - - } else { - switch (ti->type) { - case (kDouble) : - if (this->Find(key)) { - *static_cast<Double_t*>(ti->var)=*(Double_t *)this->Find(key)->GetValuePointer(); - } else { - cout << "*** ERROR!!! Could not find " << key << " in the list of variables! ***" << endl; + if(this->Find(key)) { + VarType ty = this->Find(key)->GetType(); + if (ti->nelem>1) { + // it is an array, use the appropriateinterface + switch (ti->type) { + case (kDouble) : + this_cnt = GetArray(key,static_cast<Double_t*>(ti->var),ti->nelem); + break; + case (kInt) : + this_cnt = GetArray(key,static_cast<Int_t*>(ti->var),ti->nelem); + break; + default: + Error("THcParmList","Invalid type to read %s",key); + break; } - this_cnt=1; - break; - case (kInt) : - if (this->Find(key)) { - *static_cast<Int_t*>(ti->var)=*(Int_t *)this->Find(key)->GetValuePointer(); - } else { - cout << "*** ERROR!!! Could not find " << key << " in the list of variables! ***" << endl; + } else { + switch (ti->type) { + case (kDouble) : + if(ty == kInt) { + *static_cast<Double_t*>(ti->var)=*(Int_t *)this->Find(key)->GetValuePointer(); + } else if (ty == kDouble) { + *static_cast<Double_t*>(ti->var)=*(Double_t *)this->Find(key)->GetValuePointer(); + } else { + cout << "*** ERROR!!! Type Mismatch " << key << endl; + } + this_cnt=1; + + break; + case (kInt) : + if(ty == kInt) { + *static_cast<Int_t*>(ti->var)=*(Int_t *)this->Find(key)->GetValuePointer(); + } else if (ty == kDouble) { + *static_cast<Int_t*>(ti->var)=TMath::Nint(*(Double_t *)this->Find(key)->GetValuePointer()); + cout << "*** WARNING!!! Rounded " << key << " to nearest integer " << endl; + } else { + cout << "*** ERROR!!! Type Mismatch " << key << endl; + } + this_cnt=1; + break; + default: + Error("THcParmList","Invalid type to read %s",key); + break; } - this_cnt=1; - break; - default: - Error("THcParmList","Invalid type to read %s",key); - break; } + } else { + cout << "*** ERROR!!! Could not find " << key << " in the list of variables! ***" << endl; } if (this_cnt<=0) { if ( !ti->optional ) { @@ -518,11 +530,20 @@ Int_t THcParmList::ReadArray(const char* attrC, T* array, Int_t size) " which has length " << sz << endl; } if(size<sz) sz = size; + Int_t donint = 0; + if(ty == kDouble && typeid(array[0]) == typeid(Int_t)) { + donint = 1; // Use nint when putting doubles in nint + cout << "*** WARNING!!! Rounded " << attrC << " elements to nearest integer " << endl; + } for(cnt=0;cnt<sz;cnt++) { if(ty == kInt) { array[cnt] = ((Int_t*)vp)[cnt]; } else - array[cnt] = ((Double_t*)vp)[cnt]; + if(donint) { + array[cnt] = TMath::Nint(((Double_t*)vp)[cnt]); + } else { + array[cnt] = ((Double_t*)vp)[cnt]; + } } return(cnt); } diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx index b5fe6e3c851022b949553c89d03a6c36ed4428fd..f1ccda77f1890d99e1453e9b6928a2275943c675 100644 --- a/src/THcScintillatorPlane.cxx +++ b/src/THcScintillatorPlane.cxx @@ -38,7 +38,11 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name, fNegADCHits = new TClonesArray("THcSignalHit",16); fPlaneNum = planenum; fTotPlanes = planenum; - fNScinHits = 0; + fNScinHits = 0; + // + fMaxHits=53; + fpTimeHits=0; + fpTime = new Double_t [fMaxHits]; } //______________________________________________________________________________ THcScintillatorPlane::THcScintillatorPlane( const char* name, @@ -56,6 +60,11 @@ THcScintillatorPlane::THcScintillatorPlane( const char* name, fPlaneNum = planenum; fTotPlanes = totplanes; fNScinHits = 0; + // + fMaxHits=53; + fpTimeHits=0; + fpTime = new Double_t [fMaxHits]; + } //______________________________________________________________________________ @@ -66,8 +75,10 @@ THcScintillatorPlane::~THcScintillatorPlane() delete fNegTDCHits; delete fPosADCHits; delete fNegADCHits; - + delete fpTime; } + +//______________________________________________________________________________ THaAnalysisObject::EStatus THcScintillatorPlane::Init( const TDatime& date ) { // Extra initialization for scintillator plane: set up DataDest map @@ -336,18 +347,37 @@ Int_t THcScintillatorPlane::PulseHeightCorrection() ! reference particle, need to make sure this is big enough ! to accomodate difference in TOF for other particles ! Default value in case user hasn't defined something reasonable */ - Int_t i,j,index,nfound=0; - Double_t mintdc, maxtdc,tdctotime; - Double_t pos_ph[53],neg_ph[53],postime[53],negtime[53]; // the 53 should go in a param file (was hmax_scin_hits originally) - Bool_t keep_pos[53],keep_neg[53]; // are these all really needed? + Int_t i,j,index; + Double_t mintdc, maxtdc,tdctotime,toftolerance,tmin; + Double_t pos_ph[53],neg_ph[53],postime[53],negtime[53],scin_corrected_time[53]; // the 53 should go in a param file (was hmax_scin_hits originally) + Bool_t keep_pos[53],keep_neg[53],two_good_times[53]; // are these all really needed? Double_t dist_from_center,scint_center,hit_position,time_pos[100],time_neg[100],hbeta_pcent; + Int_t timehist[200],jmax,maxhit,nfound=0; // This seems as a pretty old-fashioned way of doing things. Is there a better way? + + + // protect against spam events + if (fNScinHits>1000) { + cout <<"Too many hits "<<fNScinHits<<" in this event! Skipping!\n"; + return -1; + } + // zero out histogram + for (i=0;i<200;i++) { + timehist[i]=0; + } + for (i=0;i<53;i++) { + keep_pos[i]=kFALSE; + keep_neg[i]=kFALSE; + two_good_times[i]=kFALSE; + } mintdc=((THcHodoscope *)GetParent())->GetTdcMin(); maxtdc=((THcHodoscope *)GetParent())->GetTdcMax(); + tdctotime=((THcHodoscope *)GetParent())->GetTdcToTime(); + toftolerance=((THcHodoscope *)GetParent())->GetTofTolerance(); // hbeta_pcent=(TH((THcHodoscope *)GetParent())->GetParent() // Horrible hack until I find out where to get the central beta from momentum!! GN hbeta_pcent=0.99776; - tdctotime=((THcHodoscope *)GetParent())->GetTdcToTime(); + for (i=0;i<fNScinHits;i++) { if ((((THcSignalHit*) fPosTDCHits->At(i))->GetData()>=mintdc) && (((THcSignalHit*) fPosTDCHits->At(i))->GetData()<=maxtdc) && @@ -360,7 +390,6 @@ Int_t THcScintillatorPlane::PulseHeightCorrection() postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosPhcCoeff(index)* TMath::Sqrt(TMath::Max(0.,(pos_ph[i]/((THcHodoscope *)GetParent())->GetHodoPosMinPh(index)-1))); postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosTimeOffset(index); - // cout <<postime[i]<<endl; neg_ph[i]=((THcSignalHit*) fNegADCHits->At(i))->GetData(); negtime[i]=((THcSignalHit*) fNegTDCHits->At(i))->GetData()*tdctotime; j=((THcSignalHit*)fNegTDCHits->At(i))->GetPaddleNumber(); @@ -368,25 +397,100 @@ Int_t THcScintillatorPlane::PulseHeightCorrection() negtime[i]=negtime[i]-((THcHodoscope *)GetParent())->GetHodoNegPhcCoeff(index)* TMath::Sqrt(TMath::Max(0.,(neg_ph[i]/((THcHodoscope *)GetParent())->GetHodoNegMinPh(index)-1))); negtime[i]=negtime[i]-((THcHodoscope *)GetParent())->GetHodoNegTimeOffset(index); - // cout <<"postime = "<<postime[i]<<" negtime = "<<negtime[i]<<endl; // Find hit position. If postime larger, then hit was nearer negative side. dist_from_center=0.5*(negtime[i]-postime[i])*((THcHodoscope *)GetParent())->GetHodoVelLight(index); scint_center=0.5*(fPosLeft+fPosRight); hit_position=scint_center+dist_from_center; - // cout <<fPosLeft<<" "<<fPosRight<<" hit position = "<<hit_position<<" "; hit_position=TMath::Min(hit_position,fPosLeft); - // cout<<hit_position<<" "; hit_position=TMath::Max(hit_position,fPosRight); - //cout<<hit_position<<endl; postime[i]=postime[i]-(fPosLeft-hit_position)/((THcHodoscope *)GetParent())->GetHodoVelLight(index); negtime[i]=negtime[i]-(hit_position-fPosRight)/((THcHodoscope *)GetParent())->GetHodoVelLight(index); time_pos[i]=postime[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent); time_neg[i]=negtime[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent); nfound++; - cout <<"time pos/neg = "<<time_pos[i]<<" "<<time_neg[i]<<endl; - // for (int k=0;k<200;k++) { + for (int k=0;k<200;k++) { + tmin=0.5*k; + if ((time_pos[i]> tmin) && (time_pos[i] < tmin+toftolerance)) { + timehist[k]++; + } + if ((time_neg[i]> tmin) && (time_neg[i] < tmin+toftolerance)) { + timehist[k]++; + } + } + nfound++; + } + } + // Find the bin with most hits + jmax=-1; + maxhit=0; + for (i=0;i<200;i++) { + if (timehist[i]>maxhit) { + jmax=i; + maxhit=timehist[i]; + } + } + if (jmax>=0) { + tmin=0.5*jmax; + for (i=0;i<fNScinHits;i++) { + if ((time_pos[i]>tmin) && (time_pos[i]<tmin+toftolerance)) { + keep_pos[i]=kTRUE; + } + if ((time_neg[i]>tmin) && (time_neg[i]<tmin+toftolerance)) { + keep_neg[i]=kTRUE; + } + } + } + // Resume regular tof code, now using time filer(?) from above + // Check for TWO good TDC hits + for (i=0;i<fNScinHits;i++) { + if ((((THcSignalHit*) fPosTDCHits->At(i))->GetData()>=mintdc) && + (((THcSignalHit*) fPosTDCHits->At(i))->GetData()<=maxtdc) && + (((THcSignalHit*) fNegTDCHits->At(i))->GetData()>=mintdc) && + (((THcSignalHit*) fNegTDCHits->At(i))->GetData()<=maxtdc) && + keep_pos[i] && keep_neg[i]) { + two_good_times[i]=kTRUE; + } + } // end of loop that finds tube setting time + fpTimeHits=0; + for (i=0;i<fNScinHits;i++) { + if (two_good_times[i]) { // both tubes fired + // correct time for everything except veloc. correction in order + // to find hit location from difference in TDC. + pos_ph[i]=((THcSignalHit*) fPosADCHits->At(i))->GetData(); + postime[i]=((THcSignalHit*) fPosTDCHits->At(i))->GetData()*tdctotime; + j=((THcSignalHit*)fPosTDCHits->At(i))->GetPaddleNumber(); + index=((THcHodoscope *)GetParent())->GetScinIndex(fPlaneNum-1,j); + postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosPhcCoeff(index)* + TMath::Sqrt(TMath::Max(0.,(pos_ph[i]/((THcHodoscope *)GetParent())->GetHodoPosMinPh(index)-1))); + postime[i]=postime[i]-((THcHodoscope *)GetParent())->GetHodoPosTimeOffset(index); + // + neg_ph[i]=((THcSignalHit*) fNegADCHits->At(i))->GetData(); + negtime[i]=((THcSignalHit*) fNegTDCHits->At(i))->GetData()*tdctotime; + j=((THcSignalHit*)fNegTDCHits->At(i))->GetPaddleNumber(); + index=((THcHodoscope *)GetParent())->GetScinIndex(fPlaneNum-1,j); + negtime[i]=negtime[i]-((THcHodoscope *)GetParent())->GetHodoNegPhcCoeff(index)* + TMath::Sqrt(TMath::Max(0.,(neg_ph[i]/((THcHodoscope *)GetParent())->GetHodoNegMinPh(index)-1))); + negtime[i]=negtime[i]-((THcHodoscope *)GetParent())->GetHodoNegTimeOffset(index); + // find hit position. If postime larger, then hit was nearer negative side + dist_from_center=0.5*(negtime[i]-postime[i])*((THcHodoscope *)GetParent())->GetHodoVelLight(index); + scint_center=0.5*(fPosLeft+fPosRight); + hit_position=scint_center+dist_from_center; + hit_position=TMath::Min(hit_position,fPosLeft); + hit_position=TMath::Max(hit_position,fPosRight); + postime[i]=postime[i]-(fPosLeft-hit_position)/((THcHodoscope *)GetParent())->GetHodoVelLight(index); + negtime[i]=negtime[i]-(hit_position-fPosRight)/((THcHodoscope *)GetParent())->GetHodoVelLight(index); + scin_corrected_time[i]=0.5*(postime[i]+negtime[i]); + fpTime[i]=scin_corrected_time[i]-(fZpos+(j%2)*fDzpos)/(29.979*hbeta_pcent); + /// cout <<"fptime ["<<i<<"]= "<<fpTime[i]<<endl; + fpTimeHits++; + } + else { // only one tube fired + scin_corrected_time[i]=0.0; // not a very good "flag" but there is the logical two_good_hits... + fpTime[i]=0.0; // NOTE: we're not incrementing fpTimeHits. In principle these two lines might be deleted } } + // Start time calculation, assume xp=yp=0 radians. project all + // time values to focal plane. use average for start time. return 0; } //_____________________________________________________________________________ diff --git a/src/THcScintillatorPlane.h b/src/THcScintillatorPlane.h index e3c5f4d5ca4e4d16972f396ed1f50c4ccf92f24f..f47b1ebdd721d7c954a84663b42abb492e349230 100644 --- a/src/THcScintillatorPlane.h +++ b/src/THcScintillatorPlane.h @@ -53,6 +53,8 @@ class THcScintillatorPlane : public THaSubDetector { Double_t GetPosRight(); Double_t GetPosOffset(); Double_t GetPosCenter(Int_t PaddleNo); // here we're counting from zero! + Double_t GetFpTime(Int_t index) { return fpTime[index];}; + Int_t GetFpTimeHits() { return fpTimeHits;}; TClonesArray* fParentHitList; @@ -68,6 +70,7 @@ class THcScintillatorPlane : public THaSubDetector { Int_t fNelem; /* Need since we don't inherit from detector base class */ Int_t fNScinHits; /* Number of hits in this plane */ + Int_t fMaxHits; /* maximum number of hits to be considered - useful for dimensioning arrays */ Double_t fSpacing; /* paddle spacing */ Double_t fSize; /* paddle size */ Double_t fZpos; /* z position */ @@ -98,7 +101,11 @@ class THcScintillatorPlane : public THaSubDetector { Double_t *fNegPed; Double_t *fNegSig; Double_t *fNegThresh; - + + // + Int_t fpTimeHits; // number of entries in *fpTime + Double_t *fpTime; // array with fpTimes for this scintillator plane + virtual Int_t ReadDatabase( const TDatime& date ); virtual Int_t DefineVariables( EMode mode = kDefine ); virtual void InitializePedestals( );