|
//____________________________________________________________________ // // Delta Tdc calibration module // * reference tube: // Left 37 by default // Right 31 by default // // select hits in reference tubes according to // | (adc - ped) - adc MIP) | < fAdcSel * adc MIP // // Fill histograms (TDC*TDCGain)i - (TDC*TDCGain)ref for hits i such // that adc/adcMIP > threshold // // This module can also deal with the INEL counters that were used for // the pp running. Since they only have a TDC signal, but otherwise are // like the BB (having a left and Right modules) same data structure this is the // best palce to insert the stuff. It is concievabel ADC signals may be added // for subsequent run periods. // The ref tubes here is 13. Only tubes 13-16 can be used in the present algortihm. // Due to halo background the difference spectra has often two peaks. The alogorithm // selects a fir range. The offsets are kind of hard-wired into the setup. A better way // would be to find the first 'real' peak. For ring larger than the reference tube ring // it will be the top peak, and for rings lower it will be the lower peak. // // //____________________________________________________________________ // // $Id: BrBbDeltaTdcCalModule.cxx,v 1.7 2002/08/15 16:17:45 videbaek Exp $ // $Author: videbaek $ // $Date: 2002/08/15 16:17:45 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> // #ifndef BRAT_BrBbDeltaTdcCalModule #include "BrBbDeltaTdcCalModule.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #ifndef WIN32 #include <iostream> #include <iomanip> #include <fstream> #else #include <iostream.h> #include <iomanip.h> #include <fstream.h> #endif #ifndef BRAT_BrDataTable #include "BrDataTable.h" #endif #ifndef ROOT_TString #include "TString.h" #endif #ifndef BRAT_BrBbDig #include "BrBbDig.h" #endif #ifndef BRAT_BrRunInfoManager #include "BrRunInfoManager.h" #endif #ifndef ROOT_TH1 #include "TH1.h" #endif #ifndef ROOT_TH2 #include "TH2.h" #endif #ifndef ROOT_TNtuple #include "TNtuple.h" #endif #ifndef BRAT_BrEventNode #include "BrEventNode.h" #endif #ifndef BRAT_BrUnits #include "BrUnits.h" #endif #ifndef ROOT_TF1 #include "TF1.h" #endif #if !defined ROOT_TSpectrum #include "TSpectrum.h" #endif //____________________________________________________________________ ClassImp(BrBbDeltaTdcCalModule); //____________________________________________________________________ BrBbDeltaTdcCalModule::BrBbDeltaTdcCalModule() : BrBbCalModule() { // Default constructor. DO NOT USE SetState(kSetup); SetDefaultParameters(); } //____________________________________________________________________ BrBbDeltaTdcCalModule::BrBbDeltaTdcCalModule(const Char_t* name, const Char_t* title) : BrBbCalModule(name, title) { // Named Constructor SetState(kSetup); SetDefaultParameters(); } //____________________________________________________________________ void BrBbDeltaTdcCalModule::SetDefaultParameters() { SetMinTdc(); // default is 10 ch SetMaxTdc(); // default is 3500 ch SetAdcSel(); // default is 0.1*adcGain0 SetEnergyThreshold(); // default is 0.8 SetFitWindow(); // default is 5 ns SetNtuple(); // default is false SetUseSlewing(); // default is false } //____________________________________________________________________ void BrBbDeltaTdcCalModule::DefineHistograms() { // Define histograms. They are: // <fill in here> if (GetState() != kInit) { Stop("DefineHistograms", "Must be called after Init"); return; } if (fLoadAscii || fCommitAscii) { if (Verbose() > 5) Warning("DefineHistograms", "No need to declare histos " "since we only load calibration from ascii file"); return; } TDirectory* saveDir = gDirectory; // list of histograms fHistDir = gDirectory->mkdir(Form("%s_DeltaTdc",GetName())); fHistDir->cd(); // Make histograms here fHistDir->mkdir("Tubes"); // tube histos fHistDir->cd("Tubes"); // big tubes fDtime = new TH1F* [fNoTubes]; Float_t timeLow = -10; Float_t timeHigh= -10; Int_t nbin = 800; // For the Inel calibration the values are not centered around zero // but goes from ~0, 2, 12 to 20 nsec. Comes from the position of the Inel // counters. // Thus use this asymmetric setup for dT spectra in that case. // if(!fIsBb){ timeLow = -25.; timeHigh = +35; nbin=240; } for (Int_t i = 0; i < fNoTubes; i++) { Int_t tube = i + 1; if (tube == fRef) continue; fDtime[i] = new TH1F(Form("dtime%02d_%02d", tube, fRef), Form("%s: Delta Time %d - %d", GetName(), tube, fRef), nbin, timeLow, timeHigh); // binning is range ns * 40 chan/ns } fHistDir->cd(); fhDeltaTime = new TH2F("deltaTime", Form("%s: Delta Time vs Tube",GetName()), fNoTubes, 0.5,0.5+fNoTubes, 200, timeLow, timeHigh); if (fNtuple){ if(fIsBb) fDeltaTdc = new TNtuple("bbDTime", "BB data for DeltaTdc", "tube:refTube:adc:refAdc:time:refTime:" "ag0:refAg0:fastTime:fastTube"); else fDeltaTdc = new TNtuple("InelDTime", "Inel data for DeltaTdc", "tube:refTube:time:refTime"); } // Summaries fSumDt = new TH1F("deltaTdc", Form("%s: DeltaTdc Dt", GetName()), fNoTubes, 0.5, fNoTubes + 0.5); fSumDt->SetMarkerStyle(22); fSumDt->SetMarkerSize(1.1); fSumDt->SetMarkerColor(2); gDirectory = saveDir; } //____________________________________________________________________ void BrBbDeltaTdcCalModule::Init() { // Job-level initialisation SetState(kInit); // init base class (BrBbCalibration and BrDetectorParamsBB) BrBbCalModule::Init(); BrCalibration::EAccessMode mode = BrCalibration::kRead; if (fLoadAscii) { mode = BrCalibration::kTransitional; fCalibration->Use("deltaTdc", mode, fNoTubes); } else { if (!fCalibration->GetAccessMode("tdcGain")) { mode = BrCalibration::kRead; if (DebugLevel() > 3) Warning("Init", "TDC Gains will be read from the SQL DB"); fCalibration->Use("tdcGain", mode, fNoTubes); } // if we want to save them (ascii or DB) if (fSaveAscii || fCommitAscii) { mode = BrCalibration::kWrite; fCalibration->Use("deltaTdc", mode, fNoTubes); } fValidTube = new Bool_t [fNoTubes]; for (Int_t i = 1; i <= fNoTubes; i++) fValidTube[i-1] = kTRUE; if(fIsBb) { if (!fCalibration->GetAccessMode("pedestal")) { mode = BrCalibration::kRead; if (DebugLevel() > 3) Warning("Init", "Pedestals will be read from the SQL DB"); fCalibration->Use("pedestal", mode, fNoTubes); fCalibration->Use("pedestalWidth", mode, fNoTubes); } if (!fCalibration->GetAccessMode("adcGain0")) { mode = BrCalibration::kRead; if (DebugLevel() > 3) Warning("Init", "ADC Gain 0 will be read from the SQL DB"); fCalibration->Use("adcGain0", mode, fNoTubes); } if (!fCalibration->GetAccessMode("adcGap")) { mode = BrCalibration::kRead; if (DebugLevel() > 3) Warning("Init", "ADC Gap 0 will be read from the SQL DB"); fCalibration->Use("adcGap", mode, fNoTubes); } if (!fCalibration->GetAccessMode("adcGapStart")) { mode = BrCalibration::kRead; if (DebugLevel() > 3) Warning("Init", "ADC Gap Start will be read from the SQL DB"); fCalibration->Use("adcGapStart", mode, fNoTubes); } // use slewing if required (2nd pass of delta tdc) if (fUseSlewing) { if (!fCalibration->GetAccessMode("slewK") || !fCalibration->GetAccessMode("slewP")) { mode = BrCalibration::kRead; if (DebugLevel() > 3) Warning("Init", "Slewing correction will be read from the SQL DB"); fCalibration->Use("slewK", mode, fNoTubes); fCalibration->Use("slewP", mode, fNoTubes); } } } } } //____________________________________________________________________ void BrBbDeltaTdcCalModule::Begin() { // Run-level initialisation SetState(kBegin); // if loading cal from ascii file if (fLoadAscii) { ReadAscii(); return; } // if reading calibration parameters from ascii file to commit them, // don't need to check if histos are booked if (fCommitAscii) return; if (!HistBooked()) { Abort("Begin", "Histograms must be booked for the" " calibration to work"); return; } if (!fCalibration->RevisionExists("*")) { Abort("Begin", " A calibration revision is missing, aborting", GetName()); return; } for (Int_t i = 1; i <= fNoTubes; i++) { fCalibration->SetDeltaTdc(i, BrBbCalibration::kCalException); if(fIsBb){ Float_t tdg = fCalibration->GetTdcGain(i); Float_t ped = fCalibration->GetPedestal(i); Float_t ag0 = fCalibration->GetAdcGain0(i); Float_t gap = fCalibration->GetAdcGap(i); Float_t gps = fCalibration->GetAdcGapStart(i); // if bad cal if (!IsValid(ped) || !IsValid(ag0) || !IsValid(tdg) || !IsValid(gap) || !IsValid(gps)) { fValidTube[i - 1] = kFALSE; if (Verbose() > 2) cout << GetName() << " *** Tube " << i << " won't be used (bad ped or adc calibration)" << endl; continue; } } else { Float_t tdg = fCalibration->GetTdcGain(i); // if bad cal if (!IsValid(tdg)) { fValidTube[i - 1] = kFALSE; if (Verbose() > 2) cout << GetName() << " *** Tube " << i << " won't be used (bad ped or adc calibration)" << endl; continue; }; } } } //____________________________________________________________________ void BrBbDeltaTdcCalModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // select hits in reference tubes according to // | (adc - ped) - adc MIP) | < 0.1 * adc MIP // // Fill histograms (TDC*TDCGain)i - (TDC*TDCGain)ref for hits i such // that adc/adcMIP > threshold SetState(kEvent); if (fLoadAscii || fCommitAscii) return; BrDataTable* hits = (BrDataTable*)inNode->GetDataTable(fDigName); if (!hits) { Warning("Event", "No data table for %s", GetName()); return; } if (!hits->GetEntries()) { Warning("Event", "No hits in datatable for %s", GetName()); return; } BrBbDig* hRef = FindReference(hits); // reference tube if (hRef) Process(hits, hRef); } //____________________________________________________________________ void BrBbDeltaTdcCalModule::Process(BrDataTable* hits, BrBbDig* ref) { // private method // ---------------------------------------------------- // fill profiles with time diff hit - ref vs energy // for hits above threshold // ---------------------------------------------------- // check the fastest hit // Only used for Ntuples... // Thus if not on this whole section can be skipped. // Float_t tfast = 99999; Float_t fast = 0; if (fIsBb && HistOn() && fNtuple) { for(Int_t i = 0; i < fNoTubes; i++) { BrBbDig* hit = (BrBbDig*)hits->At(i); Int_t tube = hit->GetTubeNo(); // get calibration Float_t tdg = fCalibration->GetTdcGain(tube); Float_t ped = fCalibration->GetPedestal(tube); Float_t ag0 = fCalibration->GetAdcGain0(tube); Float_t gap = fCalibration->GetAdcGap(tube); Float_t gps = fCalibration->GetAdcGapStart(tube); // apply needed calibration Float_t adc = hit->GetAdc() - ped; Float_t tdc = hit->GetTdc(); if (tube >= fBigT && adc > gps) adc -= gap; if (adc/ag0 < fEnergyThreshold || tdc < fMinTdc || tdc > fMaxTdc) continue; tdc *= tdg; // check the fastest hit if (tdc < tfast) { tfast = tdc; fast = tube; } } } // ------------------------------------------ // fill Calibration histograms // Float_t ped = 0; Float_t ag0 = 0; Float_t tdg = 0; Float_t gap = 0; Float_t gps = 0; for(Int_t i = 0; i < fNoTubes; i++) { BrBbDig* hit = (BrBbDig*)hits->At(i); Int_t tube = hit->GetTubeNo(); if (!fValidTube[tube - 1]) continue; if (tube == ref->GetTubeNo()) continue; Float_t adc, tdc; Float_t slk = 0; Float_t slp = 0; // get calibration if(fIsBb) { ped = fCalibration->GetPedestal(tube); ag0 = fCalibration->GetAdcGain0(tube); tdg = fCalibration->GetTdcGain(tube); gap = fCalibration->GetAdcGap(tube); gps = fCalibration->GetAdcGapStart(tube); if (fUseSlewing) { slk = fCalibration->GetSlewK(tube); slp = fCalibration->GetSlewP(tube); } // if bad cal if (!IsValid(ped) || !IsValid(ag0) || !IsValid(tdg) || !IsValid(gap) || !IsValid(gps)) continue; // apply needed calibration adc = hit->GetAdc() - ped; tdc = hit->GetTdc(); if (tube >= fBigT && adc > gps) adc -= gap; if (adc/ag0 < fEnergyThreshold || tdc < fMinTdc || tdc > fMaxTdc) continue; } else { Float_t tdg = fCalibration->GetTdcGain(tube); if (!IsValid(tdg)) continue; // apply needed calibration tdc = hit->GetTdc(); if (tdc < fMinTdc || tdc > fMaxTdc) continue; } Float_t refTime = ref->GetTdc()*fCalibration->GetTdcGain(fRef); Float_t time; if(fIsBb) time = (tdc*fCalibration->GetTdcGain(tube) - slk/TMath::Sqrt(adc/ag0) - slp/(adc/ag0)); else time = tdc*fCalibration->GetTdcGain(tube); fDtime[tube - 1]->Fill(time - refTime); fhDeltaTime->Fill(tube, time-refTime); if(fIsBb && HistOn()) { Float_t rag0 = fCalibration->GetAdcGain0(fRef); Float_t radc = ref->GetAdc(); if (fNtuple) fDeltaTdc->Fill(tube, fRef, adc, radc, time, refTime, ag0, rag0, tfast, fast); } else if(HistOn()){ if (fNtuple) fDeltaTdc->Fill(tube, fRef , time, refTime); } } } //____________________________________________________________________ BrBbDig* BrBbDeltaTdcCalModule::FindReference(BrDataTable* hits) { // private method // ----------------------------- // pick up the reference tube hit // ----------------------------- BrBbDig* ref = 0; for(Int_t i = 0; i < fNoTubes; i++) { BrBbDig* hit = (BrBbDig*)hits->At(i); Int_t tube = hit->GetTubeNo(); if (tube != fRef) continue; if(fIsBb){ // get calibration Float_t ped = fCalibration->GetPedestal(tube); Float_t ag0 = fCalibration->GetAdcGain0(tube); Float_t tdg = fCalibration->GetTdcGain(tube); if (!IsValid(ped) || !IsValid(ag0)) { Abort("FindReference", "Bad calibration for reference tube %d", fRef); return 0; } Float_t adc = hit->GetAdc() - ped; Float_t tdc = hit->GetTdc(); if (TMath::Abs(adc - ag0) < fAdcSel*ag0 && tdc < fMaxTdc && tdc > fMinTdc) ref = hit; } else { // get calibration Float_t tdg = fCalibration->GetTdcGain(tube); Float_t tdc = hit->GetTdc(); if (tdc < fMaxTdc && tdc > fMinTdc) ref = hit; } } return ref; } //____________________________________________________________________ void BrBbDeltaTdcCalModule::End() { // Run-level finalisation SetState(kEnd); } //____________________________________________________________________ void BrBbDeltaTdcCalModule::Finish() { // Job-level finalisation // // This module also is useful for the pp Inel counters // The time differences can be pretty large, and in particular for the // left rings the data were plagued by considerable amount of 'background ie. // halo from the wrong direction. The correct peak was then NOT always the // one at the maximum in the spectrum (but should usualy be). // The algortim to find peaks is to use the TSpectrum class and looking for two peaks // If there are two present it will pick the proper one by // SetState(kFinish); if (fLoadAscii) return; if (fCommitAscii) { ReadAscii(); return; } if (Verbose()) cout << " *** " << GetName() << " DeltaTdc correction: " << endl; // The delta time is relative to reference tube. fCalibration->SetDeltaTdc(fRef, 0); fSumDt->Fill(fRef, 0); TSpectrum* spec; if(!fIsBb) spec = new TSpectrum(2,2); TF1* fit = new TF1("fit", "gaus", 0., 1000.); if (Verbose()) cout << " Tube | FitMean | FitSigma | Mean | RMS |" << " DMean | DSigma " << endl; for (Int_t i = 0; i < fNoTubes; i++) { Int_t tube = i + 1; if (tube == fRef) continue; if (!fValidTube[tube-1]) { cout << " bad tube: " << tube << endl; continue; } if (fDtime[i]->GetEntries() == 0) continue; // fit tube pedestals with gaussian Float_t mean = Float_t(fDtime[i]->GetMean()); Float_t width = Float_t(fDtime[i]->GetRMS()); if(!fIsBb) { spec->Search(fDtime[i], 4.); if(spec->GetNPeaks() == 1) { Int_t bin = fDtime[i]->GetMaximumBin(); mean = fDtime[i]->GetBinCenter(bin); } else { Int_t ip=1; if(tube > fRef) ip=0; mean = spec->GetPositionX()[ip]; } fDtime[i]->GetXaxis()->SetRangeUser(mean-3., mean+3.); } fit->SetParameters(fDtime[i]->GetMaximum(), mean, width); if(fDtime[i]->GetMaximum()< 10) fDtime[i]->Fit("fit", "Q0L", "", mean - 2*width, mean + 2*width); else fDtime[i]->Fit("fit", "Q0", "", mean - 2*width, mean + 2*width); Double_t meanFit = fit->GetParameter(1); Double_t sigmaFit = fit->GetParameter(2); fSumDt->Fill(tube, meanFit); if (Verbose()) cout << setw(3) << tube << setw(10) << setprecision(4) << meanFit << setw(10) << setprecision(4) << sigmaFit << setw(10) << setprecision(4) << mean << setw(10) << setprecision(4) << width << setw(10) << setprecision(4) << meanFit - mean << setw(10) << setprecision(4) << sigmaFit - width << endl; if (TMath::Abs(meanFit - mean) > 20) { meanFit = mean; sigmaFit = width; if (Verbose()) cout << " *** Deviation too large, probably bad fit" << " Take histo MEAN and RMS instead" << endl; } fCalibration->SetDeltaTdc(tube, meanFit); } // save stuff to ascii file if required if (fSaveAscii) SaveAscii(); } //____________________________________________________________________ void BrBbDeltaTdcCalModule::SaveAscii() { // save pedestal to ascii file BrRunInfoManager* runMan = BrRunInfoManager::Instance(); Int_t* runs = runMan->GetRunNumbers(); Int_t nrun = runMan->GetNumberOfRuns(); BrBbCalModule::SaveAscii(); ofstream file(fCalibFile.Data(), ios::out); file << "****************************************** " << endl; file << "* Calibration for BB counter " << GetName() << endl; file << "* DeltaTdc calibration, reference tube: " << fRef << endl; file << "* Used events from run(s) "; for (Int_t r = 0; r < nrun; r++) file << runs[r] << " "; file << endl; file << "****************************************** " <<endl; file << "*" << endl; file << "* tube | Delta Tdc | " << endl; file << "* -------------------- " << endl << endl; for (Int_t i = 0; i < fNoTubes; i++) { Int_t tube = i + 1; file << setw(5) << tube << setw(13) << fCalibration->GetDeltaTdc(tube) << endl; } file << "* ------------------------------------" << endl << endl; cout << GetName() << " DeltaTdc parameters saved to file " << fCalibFile.Data() << endl << endl; } //____________________________________________________________________ void BrBbDeltaTdcCalModule::ReadAscii() { BrBbCalModule::ReadAscii(); ifstream file(fCalibFile.Data(), ios::in); if (!file) { Abort("ReadAscii", "File %s was not found", fCalibFile.Data()); return; } Float_t dt; Int_t tube; Char_t comment[256]; file.getline(comment, 256); while(comment[0] == '*') { file.getline(comment, 256); if (DebugLevel() > 5) cout << comment << endl; } for (Int_t i = 1; i <= fNoTubes; i++) { file >> tube >> dt; if (DebugLevel() > 5) cout << setw(4) << tube << setw(12) << dt << endl; fCalibration->SetDeltaTdc(tube, dt); } fCalibration->SetComment("deltaTdc", fComment.Data()); } //____________________________________________________________________ void BrBbDeltaTdcCalModule::Print(Option_t* option) const { // Print module information // See BrModule::Print for options. // In addition this module defines the Option: // <fill in here> TString opt(option); opt.ToLower(); BrModule::Print(option); if (opt.Contains("d")) cout << endl << " Original author: Djamel Ouerdane" << endl << " Last Modifications: " << endl << " $Author $" << endl << " $Date $" << endl << " $Revision $ " << endl << endl << "-------------------------------------------------" << endl; } //____________________________________________________________________ // // $Log: BrBbDeltaTdcCalModule.cxx,v $ // Revision 1.7 2002/08/15 16:17:45 videbaek // Improved method for Inel counters that may have multiple peaks // in the dT spectrum and even not aligned.l // // Revision 1.6 2002/05/07 15:19:04 bjornhs // Changed a histogram upper limit from -10 to +10 - it made no sense to have // upperlimit=lowerlimit in the delta time plot... // // Revision 1.5 2002/03/21 15:02:16 ouerdane // added fComment member to base class module and method SetComment so that the user can set comments about the calibration at commit time // // Revision 1.4 2002/03/20 19:30:52 videbaek // Add the INL and INR counters to be dealt with for calibration and timing // In the 2002 pp run there were not really of BB kind since they only have // td information, but it is most convinientn to deal with them as BB calibration // data. // // Revision 1.3 2001/11/05 06:56:32 ouerdane // added the possibility to use the slewing correction in the delta tdc calculation // // Revision 1.2 2001/10/23 20:47:43 ouerdane // Updated calibration modules for more efficiency, details are mainly technical // // Revision 1.1 2001/10/16 02:08:49 ouerdane // Added the delta tdc calibration. Needed for slewing correction // since we want to select within 2 sigma dt (remove outlier hits // that distort the slewing profile). // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|