|
//____________________________________________________________________ // // Gain calibration module for C1. The gain is calibrated by finding // the peak in the adc spectrum that are generated by tracks // with a momentum between 4 and 6 hitting in the middle of a tube. // // In the inNode passed to the event method there should be FFSracks // (in a BrDataTable) and raw C1 data (DigC1). // // When gain calibrating C1 with this module it is very important // that the C1 geometry is aligned properly - if it's not, the tube // might get a smaller fraction of the generated Cherenkov photons // (compared to if the track hits right in the middle of the tube), // and the gain calibration constant will be wrong. // //____________________________________________________________________ // // $Id: BrC1AdcGainCalModule.cxx,v 1.3 2001/11/05 06:57:37 ouerdane Exp $ // $Author: ouerdane $ // $Date: 2001/11/05 06:57:37 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> // #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 // ROOT Includes #ifndef ROOT_TMath #include "TMath.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #ifndef ROOT_TString #include "TString.h" #endif // Brat stuff #ifndef BRAT_BrC1AdcGainCalModule #include "BrC1AdcGainCalModule.h" #endif #ifndef BRAT_BrGeometryDbManager #include "BrGeometryDbManager.h" #endif #ifndef BRAT_BrParameterDbManager #include "BrParameterDbManager.h" #endif #ifndef BRAT_BrDataTable #include "BrDataTable.h" #endif #ifndef BRAT_BrSpectrometerTracks #include "BrSpectrometerTracks.h" #endif #ifndef BRAT_BrC1Dig #include "BrC1Dig.h" #endif #ifndef BRAT_BrRunInfoManager #include "BrRunInfoManager.h" #endif #ifndef BRAT_BrRunInfo #include "BrRunInfo.h" #endif //________________________________________________________________________ ClassImp(BrC1AdcGainCalModule); //________________________________________________________________________ BrC1AdcGainCalModule::BrC1AdcGainCalModule() : BrChkvCalModule() { // // default constructor // DO NOT USE SetState(kSetup); SetMaxGainLimit(); // default is 400 SetMinGainLimit(); // default is 20 SetMomentumCut(); // default is 4, 6 (GeV/c) SetWindowWidth(); // default is 3 (cm) } //________________________________________________________________________ BrC1AdcGainCalModule::BrC1AdcGainCalModule(Char_t *Name, Char_t *Title) : BrChkvCalModule(Name, Title) { // Named Constructor SetState(kSetup); SetMaxGainLimit(); // default is SetMinGainLimit(); // default is SetMomentumCut(); // default is 4, 6 (GeV/c) SetWindowWidth(); // default is 3 (cm) } //________________________________________________________________________ BrC1AdcGainCalModule::~BrC1AdcGainCalModule() { // // destructor // } //________________________________________________________________________ void BrC1AdcGainCalModule::DefineHistograms() { // // Book histograms (called from BrModule::Book) // 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 or commit from ascii file"); return; } if (DebugLevel() > 1) cout << "In BrC1AdcGainCalModule::DefineHistograms() " << endl; TDirectory* saveDir = gDirectory; // Store a pointer to the histogram directory. This is needed, // only because we make histograms at Finish time, and we are not // garantied that we are in the histogram file at that point. fHistDir = gDirectory->mkdir(Form("%s_AdcGain","C1")); fHistDir->mkdir("ADC"); fHistDir->cd("ADC"); // adc histos fAdc = new TH1F* [fParamsChkv->GetNoTubes()]; for(Int_t t = 1; t <= fParamsChkv->GetNoTubes(); t++) fAdc[t] = new TH1F(Form("AdcTube%02d", t), Form("%s Tube %d ADC", "C1", t), 100, 0, 1000); // adc cal histos fAdcCal = new TH1F* [fParamsChkv->GetNoTubes()]; // other histos for checking // gain vs tube number (summary histos) fHistDir->cd(); fSum = new TH1F("gain", "ADC Gain summary", fParamsChkv->GetNoTubes(), 0.5, fParamsChkv->GetNoTubes() + 0.5); fSum->SetMarkerStyle(22); fSum->SetMarkerSize(1.2); fSumW = new TH1F("width", "ADC Gain Sigma summary", fParamsChkv->GetNoTubes(), 0.5, fParamsChkv->GetNoTubes() + 0.5); fSumW->SetMarkerStyle(22); fSumW->SetMarkerSize(1.2); fChi2 = new TH1F("chi2", "Fit Chi2 summary", fParamsChkv->GetNoTubes(), 0.5, fParamsChkv->GetNoTubes() + 0.5); fChi2->SetMarkerStyle(22); fChi2->SetMarkerSize(1.2); gDirectory = saveDir; } //_______________________________________________________________________ void BrC1AdcGainCalModule::Init(){ // Job-level initialisation SetState(kInit); // // Called once per session // Get detectorParameters. Inform the parameterElementManeger about // database tables to be used (i.e. filled). if(DebugLevel() > 1) cout << "Entering Init of BrC1AdcGainCalModule" << endl; BrChkvCalModule::Init(); // fCalibration is a member of the base class and has already been // registered in BrC1CalModule::Init(); Int_t nTubes = fParamsChkv->GetNoTubes(); BrCalibration::EAccessMode mode = BrCalibration::kRead; // check if we want to load adc gain cal from ascii file if (fLoadAscii) { mode = BrCalibration::kTransitional; fCalibration->Use("adcGain", mode, nTubes); } else { if (!fCalibration->GetAccessMode("pedestal") || !fCalibration->GetAccessMode("pedestalWidth")) { mode = BrCalibration::kRead; if (DebugLevel() > 3) Warning("Init", "Pedestals will be read from the SQL database"); fCalibration->Use("pedestal", mode, nTubes); fCalibration->Use("pedestalWidth", mode, nTubes); } // if we want to save them (ascii or DB) if (fSaveAscii || fCommitAscii) { mode = BrCalibration::kWrite; fCalibration->Use("adcGain", mode, nTubes); } } BrGeometryDbManager *geoDb = BrGeometryDbManager::Instance(); if (fCommitAscii || fLoadAscii) return; if (!geoDb) { Abort("Init ", "Couldn't instantiate geometry manager"); return; } fC1Volume = (BrDetectorVolume*)geoDb->GetDetectorVolume("BrDetectorVolume","C1"); fT2Volume = (BrDetectorVolume*)geoDb->GetDetectorVolume("BrDetectorVolume","T2"); fC1BackPlane = fC1Volume->GetBackPlane(); } //________________________________________________________________________ void BrC1AdcGainCalModule::Begin(){ // // Called once per run // // Run-level initialisation SetState(kBegin); if (fLoadAscii) { ReadAscii(); return; } if (!HistBooked()) if (!fCommitAscii) { Abort("Begin", "Need histos for calibration"); return; } // check if a pedestal revision exists if (!fCalibration->RevisionExists("*")) { Abort("Begin", "Could not find a pedestal calibration revision!"); return; } for (Int_t t = 1; t <= fParamsChkv->GetNoTubes(); t++) { fValidTube[t-1] = kTRUE; fCalibration->SetAdcGain(t, BrChkvCalibration::kCalException); if (!IsValid(fCalibration->GetPedestal(t) || !IsValid(fCalibration->GetPedestalWidth(t)))) fValidTube[t-1] = kFALSE; } } //________________________________________________________________________ void BrC1AdcGainCalModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // Procedure: // 1. Get FFS track table and check if there's any tracks // 1. Get DigC1 object. // 2. Loop over tracks and select tracks in momentum range // 3. Check if the track hits in the middle of a tube // 4. Subtract pedestal and fill histogram // Per event method SetState(kEvent); if (fCommitAscii || fLoadAscii) return; if (DebugLevel() > 2) cout << "Entering Event() in BrC1AdcGainCalModule for " << "C1" << endl; // get ffs track table BrDataTable* ffsTrackTable = (BrDataTable*) inNode->GetDataTable("FfsTracks"); if (!ffsTrackTable) { ffsTrackTable = (BrDataTable*) outNode->GetDataTable("FfsTracks"); if (!ffsTrackTable) { if (Verbose()>10) Warning("Event","Could not find FFS track table."); return; } } // only look at 1 track events if (ffsTrackTable->GetEntries()!=1) return; // getting digC1 data BrC1Dig* digC1 = (BrC1Dig*) inNode->GetObject("DigC1"); if (!digC1) { digC1 = (BrC1Dig*) outNode->GetObject("DigC1"); if (!digC1) { if (DebugLevel() > 5) Warning("Event", "No DigC1 in event node"); return; } } // loop over tracks for (Int_t i=0; i<ffsTrackTable->GetEntries(); i++){ BrFfsTrack* track = (BrFfsTrack*)ffsTrackTable->At(i); // only look at tracks with p in momentum cut range if ((track->GetMomentum() < fMomentumCutHigh)&& (track->GetMomentum() > fMomentumCutLow)) { // get x position of intersection with c1 backplane Float_t hitX = PointToC1BackPlane(track->GetBackTrack()).GetX(); Float_t hitY = PointToC1BackPlane(track->GetBackTrack()).GetY(); // find tube that is hit Int_t hitTube = fParamsChkv->GetClosestTube(hitX,hitY); Float_t dX = hitX - fParamsChkv->GetTubePosition(hitTube).GetX(); Float_t dY = hitY - fParamsChkv->GetTubePosition(hitTube).GetY(); // tube should be hit right in the middle if ( (TMath::Abs(dX)<fWindowWidth) && (TMath::Abs(dY)<fWindowWidth) ) { if (DebugLevel()>10) { cout << "Track with momentum " << track->GetMomentum() << " points to the middle of tube " << hitTube << "." << endl; } //get pedestal for tube t Float_t ped = fCalibration->GetPedestal(hitTube); if (ped==BrChkvCalibration::kCalException) continue; // substract pedestal Double_t adc = digC1->GetAdc(hitTube) - ped; fAdc[hitTube]->Fill(adc); } // end if track in middle of tube } // end if track has right momentum } // end of ffs tracks loop } //____________________________________________________________________ BrVector3D BrC1AdcGainCalModule::PointToC1BackPlane(BrDetectorTrack* t2track) { // Returns the position (in local C1 coordinates) of the // t2 track intersection with C1 backplane. BrLine3D trackLineG = fT2Volume->LocalToGlobal(t2track->GetTrackLine()); BrVector3D c1HitG = fC1BackPlane.GetIntersectionWithLine(trackLineG); BrVector3D c1HitL; fC1Volume->GlobalToLocal(c1HitG,c1HitL,0); return c1HitL; } //________________________________________________________________________ void BrC1AdcGainCalModule::End(){ // // Called once at end of each run // // Run-level finalisation SetState(kEnd); } //________________________________________________________________________ void BrC1AdcGainCalModule::Finish() { // Job-level finalisation SetState(kFinish); //-------------------------------------------------- // fit ADC distribution to get the peak channel. // this channel will be the gain passed to the calibration // parameter element //-------------------------------------------------- // if load ascii mode if (fLoadAscii) return; // if commit mode if (fCommitAscii) { ReadAscii(); return; } TF1* fit = new TF1("Fit", "gaus", 0., 500.); fit->SetParNames("A","GainMean","GainWidth"); if (Verbose() > 2) cout << "C1" << " ADC histos ready " << endl << " starting fitting procedure with gaus" << endl << endl << " **** Fit parameters **** " << endl << endl << " Tube | Gain Mean | Gain sigma | Chi2 " << endl << endl; // Go to the histogram directory, as stored in the member fHistDir. TDirectory* savDir = gDirectory; fHistDir->cd("ADC"); //--- tubes for (Int_t i = 1; i <= fParamsChkv->GetNoTubes(); i++) { EvaluateGain(fAdc[i], fit, i); MakeAdcCalHist(i); } fCalibration->SetComment("adcGain", "Generated by BrC1AdcGainCalModule: " "fit with gaus"); // Go back to old directory gDirectory = savDir; if (fSaveAscii) SaveAscii(); } //____________________________________________________________________ void BrC1AdcGainCalModule::EvaluateGain(TH1F* h, TF1* fit, Int_t tube) { //-------------------------------------------------- // fit ADC distribution to get the peak channel. // this channel will be the gain passed to the calibration // parameter element //-------------------------------------------------- Double_t maxStat = h->GetMaximum(); Double_t meanStat = h->GetMean(); // cout << meanStat << " " << maxStat << endl; // note that you can choose another parameter than maxStat // set parameters fit->SetParameter(0, maxStat); fit->SetParameter(1, meanStat); fit->SetParameter(2, 30); // set limits to adc gain fit->SetParLimits(4, fMinGainLimit, fMaxGainLimit); // should maybe set a range for the fit... h->Fit(fit->GetName(), "Q0"); Double_t chi2 = fit->GetChisquare() / fit->GetNDF(); Double_t gain = fit->GetParameter("GainMean"); Double_t gainSigma = fit->GetParameter("GainWidth"); // check that result is not too crazy if (gain < fMinGainLimit || gain > fMaxGainLimit) { gain = BrChkvCalibration::kCalException; gainSigma = BrChkvCalibration::kCalException; } if (Verbose() > 1) Printf(" Tube %d %03f %03f %03f", tube , gain, gainSigma, chi2); fCalibration->SetAdcGain(tube, gain); fSum->Fill(tube, gain); fSumW->Fill(tube, gainSigma); fChi2->Fill(tube, chi2); } //_____________________________________________________________________ void BrC1AdcGainCalModule::MakeAdcCalHist(Int_t tube) { // PRIVATE METHOD: // Multiply ADC by gain (divide by ADC at peak) // and fill calibrated adc histogram Int_t minBin = fAdc[tube]->GetXaxis()->GetFirst(); Int_t maxBin = fAdc[tube]->GetXaxis()->GetLast(); Axis_t xmin = fAdc[tube]->GetXaxis()->GetXmin()/fCalibration->GetAdcGain(tube); Axis_t xmax = fAdc[tube]->GetXaxis()->GetXmax()/fCalibration->GetAdcGain(tube); Float_t factor = (xmax - xmin)/Float_t(fAdc[tube]->GetNbinsX()); fAdcCal[tube] = new TH1F(Form("AdcCalTube%02d", tube), Form("C1 Tube %d Cal ADC", tube), maxBin - minBin, xmin, xmax); for (Int_t i = minBin; i <= maxBin; i++) { Axis_t x = i*factor; fAdcCal[tube]->Fill(x, fAdc[tube]->GetBinContent(i)); } } //____________________________________________________________________ void BrC1AdcGainCalModule::SaveAscii() { // save pedestal to ascii file BrRunInfoManager* runMan = BrRunInfoManager::Instance(); Int_t* runs = runMan->GetRunNumbers(); Int_t nrun = runMan->GetNumberOfRuns(); BrChkvCalModule::SaveAscii(); ofstream file(fCalibFile.Data(), ios::out); file << "****************************************** " << endl; file << "* Gain Calibration for C1 detector " << endl; file << "* ADC 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 << "* tube | adc gain | " << endl; file << "* ---------------------" << endl << endl; for (Int_t i = 0; i < fParamsChkv->GetNoTubes(); i++) { Int_t tube = i + 1; file << setw(4) << tube << setw(15) << fCalibration->GetAdcGain(tube) << endl; } file << "* ---------------------" << endl << endl; } //____________________________________________________________________ void BrC1AdcGainCalModule::ReadAscii() { // read ascii calibration produced by this module BrChkvCalModule::ReadAscii(); ifstream file(fCalibFile.Data(), ios::in); if (!file) { Stop("ReadAscii", "File %s was not found", fCalibFile.Data()); return; } Float_t gain; 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 <= fParamsChkv->GetNoTubes(); i++) { file >> tube >> gain; if (DebugLevel() > 5) cout << setw(4) << tube << setw(12) << gain << endl; fCalibration->SetAdcGain(tube, gain); } fCalibration->SetComment("adcGain", "Generated by BrC1AdcGainCalModule: " "fit with gaus function"); } //________________________________________________________________________ void BrC1AdcGainCalModule::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: Claus Ekman Jorgensen" << endl << " Last Modifications: " << endl << " $Author: ouerdane $" << endl << " $Date: 2001/11/05 06:57:37 $" << endl << " $Revision: 1.3 $ " << endl << endl << "-------------------------------------------------" << endl; } //____________________________________________________________________ // // $Log: BrC1AdcGainCalModule.cxx,v $ // Revision 1.3 2001/11/05 06:57:37 ouerdane // changed FFS to Ffs // // Revision 1.2 2001/10/23 20:51:23 ouerdane // Updated modules ala tof or bb, with Set[Save,Commit,Load]Ascii, etc // // Revision 1.1 2001/09/03 18:00:43 ekman // Added BrC1AdcGainCalModule. // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|