|
//____________________________________________________________________ // // TDC gain calibration module // Use events from TDC calibration run (e.g 2544) // Fill histograms and profile, check linearity of the TDCs //____________________________________________________________________ // // $Id: BrTofTdcGainCalModule.cxx,v 1.10 2002/03/21 15:04:25 ouerdane Exp $ // $Author: ouerdane $ // $Date: 2002/03/21 15:04:25 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> // #ifndef BRAT_BrTofTdcGainCalModule #include "BrTofTdcGainCalModule.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_BrTofDig #include "BrTofDig.h" #endif #ifndef BRAT_BrRunInfoManager #include "BrRunInfoManager.h" #endif #ifndef ROOT_TH2 #include "TH2.h" #endif #ifndef ROOT_TFile #include "TFile.h" #endif //____________________________________________________________________ ClassImp(BrTofTdcGainCalModule); //____________________________________________________________________ BrTofTdcGainCalModule::BrTofTdcGainCalModule() { // Default constructor. DO NOT USE SetState(kSetup); SetTimeStep(); // default is 10 ns SetMinGain(); // default is 0.020 ns / chan SetMaxGain(); // default is 0.030 ns / chan } //____________________________________________________________________ BrTofTdcGainCalModule::BrTofTdcGainCalModule(const Char_t* name, const Char_t* title) : BrTofCalModule(name, title) { // Named Constructor SetState(kSetup); SetTimeStep(); // default is 10 ns SetMinGain(); // default is 0.020 ns / chan SetMaxGain(); // default is 0.030 ns / chan } //____________________________________________________________________ void BrTofTdcGainCalModule::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; fHistDir = gDirectory->mkdir(Form("%sTdcGain",GetName())); fHistDir->cd(); // Make histograms here fHistDir->mkdir("Top"); fHistDir->mkdir("Bot"); // top tube histos fHistDir->cd("Top"); fTTdc = new TH1F* [fParamsTof->GetNoSlats()]; for (Int_t i = 1; i <= fParamsTof->GetNoSlats(); i++) fTTdc[i-1] = new TH1F(Form("ttdc%d", i), Form("%s Tdc Top tube %d", GetName(), i), 4000, 0.5, 4000.5); // bot tube histos fHistDir->cd("Bot"); fBTdc = new TH1F* [fParamsTof->GetNoSlats()]; for (Int_t i = 1; i <= fParamsTof->GetNoSlats(); i++) fBTdc[i-1] = new TH1F(Form("btdc%d", i), Form("%s Tdc Bot tube %d", GetName(), i), 4000, 0.5, 4000.5); // top tubes fHistDir->cd("Top"); fTGain = new TH1F* [fParamsTof->GetNoSlats()]; for (Int_t i = 1; i <= fParamsTof->GetNoSlats(); i++) { fTGain[i-1] = new TH1F(Form("tgain%d", i), Form("%s Tdc Gain, top tube %d", GetName(), i), 1200, 0.5, 120.5); fTGain[i-1]->SetMarkerStyle(22); fTGain[i-1]->SetMarkerSize(1.2); fTGain[i-1]->SetMarkerColor(2); } // bot tubes fHistDir->cd("Bot"); fBGain = new TH1F* [fParamsTof->GetNoSlats()]; for (Int_t i = 1; i <= fParamsTof->GetNoSlats(); i++) { fBGain[i-1] = new TH1F(Form("bgain%d", i), Form("%s Tdc gain, bot tube %d", GetName(), i), 1200, 0.5, 120.5); fBGain[i-1]->SetMarkerStyle(22); fBGain[i-1]->SetMarkerSize(1.2); fBGain[i-1]->SetMarkerColor(2); } // summaries fHistDir->cd(); // top tubes fTSum = new TH1F("topGains", Form("%s TDC Gains Top tubes", GetName()), fParamsTof->GetNoSlats(), 0.5, fParamsTof->GetNoSlats() + 0.5); fTSum->SetMarkerStyle(22); fTSum->SetMarkerSize(1.1); fTSum->SetMarkerColor(2); // bot tubes fBSum = new TH1F("botGains", Form("%s TDC Gains Bot tubes", GetName()), fParamsTof->GetNoSlats(), 0.5, fParamsTof->GetNoSlats() + 0.5); fBSum->SetMarkerStyle(22); fBSum->SetMarkerSize(1.1); fBSum->SetMarkerColor(4); gDirectory = saveDir; } //____________________________________________________________________ void BrTofTdcGainCalModule::Init() { // Job-level initialisation SetState(kInit); // init base class (BrTofCalibration and BrDetectorParamsTof) BrTofCalModule::Init(); Int_t nSlats = fParamsTof->GetNoSlats(); BrCalibration::EAccessMode mode = BrCalibration::kRead; // check if we want to load cscint and delta delay cal // from ascii file if (fLoadAscii) { mode = BrCalibration::kTransitional; fCalibration->Use("topTdcGain", mode, nSlats ); fCalibration->Use("botTdcGain", mode, nSlats); } else if (fSaveAscii || fCommitAscii) { // check if need calibration already loaded // if not, try to get them from DB mode = BrCalibration::kWrite; fCalibration->Use("topTdcGain", mode, nSlats ); fCalibration->Use("botTdcGain", mode, nSlats); } } //____________________________________________________________________ void BrTofTdcGainCalModule::Begin() { // Run-level initialisation SetState(kBegin); if (fLoadAscii) { ReadAscii(); return; } if (!HistBooked()) if (!fCommitAscii) { Abort("Begin", "Histograms must be booked for the" " calibration to work"); return; } for (Int_t s = 1; s <= fParamsTof->GetNoSlats(); s++) { fCalibration->SetTopTdcGain(s, BrTofCalibration::kCalException); fCalibration->SetBotTdcGain(s, BrTofCalibration::kCalException); } } //____________________________________________________________________ void BrTofTdcGainCalModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // Fill histograms with TDC values for top and bot tubes // note that you should use only a calibration run (e.g. 2544) // Per event method SetState(kEvent); if (fCommitAscii || fLoadAscii) return; BrDataTable* hits = (BrDataTable*)inNode->GetDataTable(Form("DigTof %s", GetName())); if (!hits) { if (DebugLevel() > 5) Warning("Event", "No data table for %s", GetName()); return; } if (!hits->GetEntries()) { if (DebugLevel() > 5) Warning("Event", "No hits in datatable for %s", GetName()); return; } for(Int_t i = 0; i < fParamsTof->GetNoSlats(); i++) { BrTofDig* hit = (BrTofDig*)hits->At(i); Int_t slat = hit->GetSlatno(); if (hit->GetTdcUp() > 1 && hit->GetTdcUp() < 4095) fTTdc[slat-1]->Fill(hit->GetTdcUp()); if (hit->GetTdcDown() > 1 && hit->GetTdcDown() < 4095) fBTdc[slat-1]->Fill(hit->GetTdcDown()); } } //____________________________________________________________________ void BrTofTdcGainCalModule::Finish() { // Job-level finalisation SetState(kFinish); // if load ascii mode if (fLoadAscii) return; // if commit mode if (fCommitAscii) { ReadAscii(); return; } // determine maxima (sharp peaks) // between channel 0 and 4000, there are 10 peaks, the time step in // the calibrations are assumed to be 10 nsec // Top tubes for (Int_t i = 1; i <= fParamsTof->GetNoSlats(); i++) { // scan histo h, find each peak // result stored in profile for final fit FindMax(fTTdc[i-1], fTGain[i-1]); fTGain[i-1]->Fit("pol1", "Q0"); TF1* pol = (TF1*)fTGain[i-1]->GetListOfFunctions()->At(0); if (Verbose() > 3) cout << fTGain[i-1]->GetName() << " : " << pol->GetParameter(1) << " +/- " << pol->GetParError(1) << endl; if (1/pol->GetParameter(1) < fMinGain || 1/pol->GetParameter(1) > fMaxGain) fCalibration->SetTopTdcGain(i, BrTofCalibration::kCalException); else fCalibration->SetTopTdcGain(i, 1/pol->GetParameter(1)); } //---------- // Bot tubes for (Int_t i = 1; i <= fParamsTof->GetNoSlats(); i++) { // scan histo h, find each peak // result stored in profile for final fit FindMax(fBTdc[i-1], fBGain[i-1]); fBGain[i-1]->Fit("pol1", "Q0"); TF1* pol = (TF1*)fBGain[i-1]->GetListOfFunctions()->At(0); if (Verbose() > 3) cout << fBGain[i-1]->GetName() << " : " << pol->GetParameter(1) << " chan/ns +/- " << pol->GetParError(1) << endl; if (1/pol->GetParameter(1) < fMinGain || 1/pol->GetParameter(1) > fMaxGain) fCalibration->SetBotTdcGain(i, BrTofCalibration::kCalException); else fCalibration->SetBotTdcGain(i, 1/pol->GetParameter(1)); } // time to summarize: for (Int_t i = 1; i <= fParamsTof->GetNoSlats(); i++) { fTSum->Fill(i, (Stat_t)fCalibration->GetTopTdcGain(i)*1000); fBSum->Fill(i, (Stat_t)fCalibration->GetBotTdcGain(i)*1000); } // save stuff to ascii file if required if (fSaveAscii) SaveAscii(); } //____________________________________________________________________ void BrTofTdcGainCalModule::FindMax(TH1F* h, TH1F* p) { // ---------------- // scan histo h, find each peak // result stored in profile for final fit // ---------------- Int_t bin = 0; // running bin Double_t tdcAtMax = 0; // tdc channel at max Int_t noPeak = 0; while (1) { // check if we are out of bound if (bin >= 4000) break; // increment bin to peek the start of the peak while (h->GetBinContent(bin) == 0) { bin++; if (bin >= 4000) break; } // recheck if we are out of bound if (bin >= 4000) break; if (DebugLevel() > 10) cout << GetName() << " TdcGain cal: bin is: " << bin << " with content: " << h->GetBinContent(bin) << endl; // now we are at start of peak // the width is less than 50 channels // ---- get channel with max stat Int_t bRange = 50; Stat_t max = 0; for (Int_t j = bin; j < bin + bRange; j++) max = TMath::Max(max, h->GetBinContent(j)); // get tdc at max for (Int_t j = bin; j < bin + bRange; j++) if (h->GetBinContent(j) == max) { tdcAtMax = h->GetBinCenter(j); break; } // // ---- Find precise mean and RMS // Stat_t stat[4]; TH1F hh(*h); hh.GetXaxis()->SetRange(tdcAtMax-10,tdcAtMax+10); hh.GetStats(stat); Double_t mean = stat[2]/stat[0]; Double_t rms = TMath::Sqrt(stat[3]/stat[0]-mean*mean); Double_t err = rms/TMath::Sqrt(stat[0]); noPeak++; if (Verbose() > 5) cout << " Peak " << setw(3) << noPeak << " at mean chan: " << setw(7) << mean << " +-" << setw(4) << err << endl; // set mean and error Float_t time = noPeak*fTimeStep; // time step in ns Double_t range = p->GetXaxis()->GetXmax() - p->GetXaxis()->GetXmin(); Double_t binX = time*p->GetNbinsX()/range; if (Verbose() > 5) cout << " Time: " << setw(5) << time << " Xrange: " << setw(5) << range << " binX: " << (Int_t)binX << endl; p->SetBinContent((Int_t)binX, mean); p->SetBinError((Int_t)binX, err); bin += 350; // go to next peak (350 is harcoded...) } return; } //____________________________________________________________________ void BrTofTdcGainCalModule::SaveAscii() { // save pedestal to ascii file BrRunInfoManager* runMan = BrRunInfoManager::Instance(); Int_t* runs = runMan->GetRunNumbers(); Int_t nrun = runMan->GetNumberOfRuns(); BrTofCalModule::SaveAscii(); ofstream file(fCalibFile.Data(), ios::out); file << "****************************************** " << endl; file << "* Calibration for Tof detector " << GetName() << endl; file << "* TDC Gain calibration " << 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 << "* slat | top TDC gain | bot TDC gain " << endl; file << "* ------------------------------------" << endl << endl; for (Int_t i = 0; i < fParamsTof->GetNoSlats(); i++) { Int_t slat = i + 1; file << setw(5) << slat << setw(13) << fCalibration->GetTopTdcGain(slat) << setw(13) << fCalibration->GetBotTdcGain(slat) << setw(13) << endl; } file << "* ------------------------------------" << endl << endl; cout << GetName() << " TDC gains saved to file " << fCalibFile.Data() << endl << endl; } //____________________________________________________________________ void BrTofTdcGainCalModule::ReadAscii() { BrTofCalModule::ReadAscii(); ifstream file(fCalibFile.Data(), ios::in); if (!file) { Stop("ReadAscii", "File %s was not found", fCalibFile.Data()); return; } Float_t tgain, bgain; Int_t slat; 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 <= fParamsTof->GetNoSlats(); i++) { file >> slat >> tgain >> bgain; if (DebugLevel() > 5) cout << setw(4) << slat << setw(12) << tgain << setw(12) << bgain << endl; fCalibration->SetTopTdcGain(slat, tgain); fCalibration->SetBotTdcGain(slat, bgain); } fCalibration->SetComment("topTdcGain", fComment.Data()); fCalibration->SetComment("botTdcGain", fComment.Data()); } //____________________________________________________________________ void BrTofTdcGainCalModule::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: ouerdane $" << endl << " $Date: 2002/03/21 15:04:25 $" << endl << " $Revision: 1.10 $ " << endl << endl << "-------------------------------------------------" << endl; } //____________________________________________________________________ // // $Log: BrTofTdcGainCalModule.cxx,v $ // Revision 1.10 2002/03/21 15:04:25 ouerdane // added fComment member to base class module and method SetComment so that the user can set comments about the calibration at commit time. Removed mean momentum stuff in slewing cal module // // Revision 1.9 2001/11/05 06:59:44 ouerdane // changed to deal with new track classes and fixed some bugs // // Revision 1.8 2001/10/08 11:28:19 cholm // Changed to use new DB access classes // // Revision 1.7 2001/10/02 01:53:28 ouerdane // Added SetSaveAscii, SetLoadAscii, SetCommitAscii and updated the way parameters are tagged for Use // // Revision 1.6 2001/07/31 09:22:10 ouerdane // Removed all references to a TList member, replaced // by histogram members, declare ntuple if fNtuple is true (set // via SetNtuple) // // Revision 1.5 2001/07/18 16:14:38 ouerdane // Updated this module: // - in Event: fill histograms only when tdc (top or bot) is > 1 and < 4095 // - added methods SetMinGain and SetMaxGain (with corresponding private members) // to set some reasonable limits to what the gains should be) // if the calculated gain is beyond this range, the // BrTofCalibration::kCalException is assigned. The range by default is // [0.023, 0.027] ns. // - added SetTimeStep and fTimeStep so that the user can set the time step of // the pulser (10 ns by default) // - made it compatible with brat2 and tested it (script available in // brahms_app/do_app/tof/mm_scripts/tofTdcGain.C and commitTdcGain.C) // // Revision 1.4 2001/07/13 20:50:16 videbaek // Improved method for defining constant. Calculate the mean and RMS for each 10 nsec peak // as Dana has been doing. // // Revision 1.3 2001/07/13 15:32:40 videbaek // The top histogram was used for the Bottom TDC gain Calculation. // The name of the Bottom histogram was wrong in fitting. resulting in bottom // gains being equel to top gains. // // Revision 1.2 2001/06/22 17:39:34 cholm // Changes t odo with data classes renaming. // // Revision 1.1.1.1 2001/06/21 14:55:06 hagel // Initial revision of brat2 // // Revision 1.1 2001/06/19 12:46:49 ouerdane // Added calibration classes and reconstruction classes. // Brat compiles but these classes haven't been tested in this // context. Be careful if you use them before I (DO) check if // all is ok. // // Note: some classes are still not included (BrTofSlewingModule, // BrTofCscintCalModule, BrTofTdcOffsetModule). Will do that after // Brat2 is available // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|