|
//____________________________________________________________________ // // BrTofSlewingCalModule: module for Tof slewing correction // // - Select tofhits and tracks that are are useful for determining // slewing. This may involve some of the following cuts // // a) Select hits matcing tracks (X and Y) // b) Cut in energy (from 0.9 to 2 MIP energy) // c) Make sure we select only pions (given by C1 or Rich information) // // - Make the selection of C1, RICH an option (see below) // - Added option to use the timing relative to the TD1 start counter for // pp running. Due to the way the slew1 (timing offset ) is defined the // difference between AA running and pp running is ~6 nsec being the flight // time between the nominal IP and the TD1 counters. // - Module as well as meaning of constants complete re-written following // the methods developed for TOFW by JH and Dana. See brief description in Event() // // // //____________________________________________________________________ // // $Id: BrTofSlewingCalModule.cxx,v 1.16 2002/08/06 19:04:30 hagel Exp $ // $Author: hagel $ // $Date: 2002/08/06 19:04:30 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> //____________________________________________________________________ // #ifndef WIN32 #include <iostream> #include <iomanip> #include <fstream> #else #include <iostream.h> #include <iomanip.h> #include <fstream.h> #endif #ifndef BRAT_BrTofSlewingCalModule #include "BrTofSlewingCalModule.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #if !defined ROOT_TObjArray #include "TObjArray.h" #endif #ifndef ROOT_TH1 #include "TH1.h" #endif #ifndef ROOT_TH2 #include "TH2.h" #endif #ifndef ROOT_TProfile #include "TProfile.h" #endif #ifndef ROOT_TF1 #include "TF1.h" #endif #ifndef ROOT_TString #include "TString.h" #endif #if !defined ROOT_TVirtualFitter #include "TVirtualFitter.h" #endif //#if !defined ROOT_TMinuit //#include "TMinuit.h" //#endif #ifndef BRAT_BrEventNode #include "BrEventNode.h" #endif #ifndef BRAT_BrDataTable #include "BrDataTable.h" #endif #ifndef BRAT_BrTofDig #include "BrTofDig.h" #endif #ifndef BRAT_BrBbVertex #include "BrBbVertex.h" #endif #ifndef BRAT_BrBbRdo #include "BrBbRdo.h" #endif #ifndef BRAT_BrRunInfoManager #include "BrRunInfoManager.h" #endif #ifndef BRAT_BrUnits #include "BrUnits.h" #endif #ifndef BRAT_BrCalibration #include "BrCalibration.h" #endif #ifndef BRAT_BrTofTrackMatch #include "BrTofTrackMatch.h" #endif #ifndef BRAT_BrTd1TrackMatch #include "BrTd1TrackMatch.h" #endif #ifndef BRAT_BrTMrsFTrackMatch #include "BrTMrsFTrackMatch.h" #endif #ifndef BRAT_BrSpectrometerTracks #include "BrSpectrometerTracks.h" #endif #ifndef BRAT_BrPid #include "BrPid.h" #endif #ifndef BRAT_BrC1Pid #include "BrC1Pid.h" #endif #ifndef BRAT_BrRichPid #include "BrRichPid.h" #endif #ifndef BRAT_BrTd1Rdo #include "BrTd1Rdo.h" #endif #include <vector> // Use vector instead of the constructs like // BtTofDig* tofhit[nmch] ; that is not standard C++ but happens to work // with gcc // // Set values of static members //____________________________________________________________________ BrTofSlewingCalModule* BrTofSlewingCalModule::fgInstance = 0; //____________________________________________________________________ //____________________________________________________________________ ClassImp(BrTofSlewingCalModule); //____________________________________________________________________ BrTofSlewingCalModule::BrTofSlewingCalModule() : BrTofCalModule() { // Default constructor. DO NOT USE SetState(kSetup); SetDefaultParameters(); fSlewHitsTable=0; } //____________________________________________________________________ BrTofSlewingCalModule::BrTofSlewingCalModule(const Char_t* name, const Char_t* title) : BrTofCalModule(name, title) { // Named Constructor fSlewHitsTable=0; SetState(kSetup); SetDefaultParameters(); } //____________________________________________________________________ BrTofSlewingCalModule::~BrTofSlewingCalModule (){ if(fSlewHitsTable) delete [] fSlewHitsTable; } //____________________________________________________________________ void BrTofSlewingCalModule::SetDefaultParameters() { // Set default params for slewing correction SetMaxTofDE(); // default is 2 MIPs SetC1Threshold(); // default is 1 SetRichPionMass(); // default is theoretical pion mass SetRichMassCut(); // default is 0.05 GeV/c^2 SetUseChkv(); // default require cherenkovs SetUseMinuit(); // SetTimeResolution(); // SetVaryDt(); SetMomentumCut(); SetUseBbVertex(); SetUseTd1Time(); // default is false (pp only) SetUseTMrsTime(); // default is false (pp only) SetTimeCut(); fNominalLength = 433.5; fFitMode = 0; } //____________________________________________________________________ void BrTofSlewingCalModule::DefineHistograms() { // Define histograms. They are: // <fill in here> // The slewing histograms are for the 'old' methods done in ADC(cal) vs time // and for the (new) minuit fitting done with the ADC(raw) vs time. // 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; } const Int_t nSlats = fParamsTof->GetNoSlats(); TDirectory* saveDir = gDirectory; fHistDir = gDirectory->mkdir(Form("%s_Slewing",GetName())); fHistDir->cd(); // Make histograms here fHistDir->mkdir("Bot"); fHistDir->mkdir("Top"); // histos for calib: Float_t tmin = -3 ; Float_t tmax = +3.; Int_t bins = TMath::Nint((tmax - tmin)/fTimeResolution); // 0.085 estimated tof res. Float_t el = 0.5; Float_t eh = 2.5; if(fUseMinuit){ el = 200; eh = 3200; } fHistDir->cd("Top"); fTSlew = new TProfile* [nSlats]; fTOff = new TH2F* [nSlats]; fTOffCorr = new TH2F* [nSlats]; for (Int_t i = 1; i <= nSlats; i++) { if(!fMinuit) fTSlew[i-1] = new TProfile(Form("top%03d", i), "", 50, el, eh, tmin, tmax); fTOff[i-1] = new TH2F(Form("ttime%03d", i), "", 100, el, eh, bins, tmin, tmax); fTOffCorr[i-1] = new TH2F(Form("ttimeCorr%03d", i), "", 100, el, eh, bins, tmin, tmax); } fHistDir->cd("Bot"); fBSlew = new TProfile* [nSlats]; fBOff = new TH2F* [nSlats]; fBOffCorr = new TH2F* [nSlats]; for (Int_t i = 1; i <= nSlats; i++) { if(!fMinuit) fBSlew[i-1] = new TProfile(Form("bot%03d", i), "", 50, el, eh, tmin, tmax); fBOff[i-1] = new TH2F(Form("btime%03d", i), "", 100, el, eh, bins, tmin, tmax); fBOffCorr[i-1] = new TH2F(Form("btimeCorr%03d", i), "", 100, el, eh, bins, tmin, tmax); } fHistDir->cd(); fTSlew1 = new TH1F("tslew1", Form("%s top tubes: slewing 1", GetName()), nSlats, 0.5, nSlats + 0.5); fTSlew2 = new TH1F("tslew2", Form("%s top tubes: slewing 2", GetName()), nSlats, 0.5, nSlats + 0.5); fBSlew1 = new TH1F("bslew1", Form("%s top tubes: slewing 1", GetName()), nSlats, 0.5, nSlats + 0.5); fBSlew2 = new TH1F("bslew2", Form("%s top tubes: slewing 2", GetName()), nSlats, 0.5, nSlats + 0.5); fTDiff = new TH2F("timediffs", Form("%s All tubes: timediffs", GetName()), nSlats, 0.5, nSlats + 0.5, 100, -2, 2); fTDiffCorr = new TH2F("timediffsCorr", Form("%s All tubes: timediffs", GetName()), nSlats, 0.5, nSlats + 0.5, 100, -2, 2); gDirectory = saveDir; } //____________________________________________________________________ void BrTofSlewingCalModule::Init() { // Job-level initialisation SetState(kInit); //--------------- // base class initialization (register calibration parameters) 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("topSlewPar1", mode, nSlats); fCalibration->Use("topSlewPar2", mode, nSlats); fCalibration->Use("botSlewPar1", mode, nSlats); fCalibration->Use("botSlewPar2", mode, nSlats); } else { // check if need calibration already loaded // if not, try to get them from DB if (!fCalibration->GetAccessMode("topPedestal") || !fCalibration->GetAccessMode("botPedestal")) { if (DebugLevel() > 3) Info(3,"Init", "Pedestals will be read from DB"); mode = BrCalibration::kRead; fCalibration->Use("topPedestal", mode, nSlats); fCalibration->Use("botPedestal", mode, nSlats); } if (!fCalibration->GetAccessMode("topAdcGain") || !fCalibration->GetAccessMode("botAdcGain")) { if (DebugLevel() > 3) Info(3,"Init", "Adc Gains will be read from DB"); mode = BrCalibration::kRead; fCalibration->Use("topAdcGain", mode, nSlats); fCalibration->Use("botAdcGain", mode, nSlats); } if (!fCalibration->GetAccessMode("topTdcGain") || !fCalibration->GetAccessMode("botTdcGain")) { if (DebugLevel() > 3) Info(3,"Init", "Tdc Gains will be read from DB"); mode = BrCalibration::kRead; fCalibration->Use("topTdcGain", mode, nSlats); fCalibration->Use("botTdcGain", mode, nSlats); } if (!fCalibration->GetAccessMode("timeOffset")) Info(3,"Init", "TimeOffset be read from DB");{ mode = BrCalibration::kRead; fCalibration->Use("timeOffset", mode, nSlats); } /* if (!fCalibration->GetAccessMode("effSpeedOfLight")) { Info(3,"Init", "Cscint will be read from DB"); mode = BrCalibration::kRead; fCalibration->Use("effSpeedOfLight", mode, nSlats); } */ // if we want to save calibration (ascii or DB) if (fSaveAscii || fCommitAscii) { mode = BrCalibration::kWrite; fCalibration->Use("topSlewPar1", mode, nSlats); fCalibration->Use("topSlewPar2", mode, nSlats); fCalibration->Use("botSlewPar1", mode, nSlats); fCalibration->Use("botSlewPar2", mode, nSlats); } } // geometry if (!fCommitAscii && !fLoadAscii) BrTofCalModule::InitGeo(); // Setup the proper Nominal length if(fInFfs){ if(fUseTd1Time) fNominalLength=680; else fNominalLength = 900; } if(fInBfs){ if(fUseTd1Time) fNominalLength=1580; else fNominalLength = 1800; } if(fInMrs) { if(fUseTMrsTime) fNominalLength = 380.; else fNominalLength = 433.5; } // Initialize fitter if Minuit is to be used. if(fUseMinuit){ fMinuit = TVirtualFitter::Fitter(this, 5); // TMinuit(5); fSlewHitsTable = new TObjArray [nSlats]; } } //____________________________________________________________________ void BrTofSlewingCalModule::Begin() { // Run-level initialisation SetState(kBegin); if (fLoadAscii) { ReadAscii(); return; } if (!HistBooked()) if (!fCommitAscii) { Abort("Begin", "Need histos for analysis"); return; } // check if we got parameter revisions: if (!fCalibration->RevisionExists("*")) { Abort("Begin", "Some calibration revision are missing!"); return; } for (Int_t s = 1; s <= fParamsTof->GetNoSlats(); s++) { fCalibration->SetTopSlewPar1(s, BrTofCalibration::kCalException); fCalibration->SetTopSlewPar2(s, BrTofCalibration::kCalException); fCalibration->SetBotSlewPar1(s, BrTofCalibration::kCalException); fCalibration->SetBotSlewPar2(s, BrTofCalibration::kCalException); fValidSlat[s-1] = kTRUE; if (!IsValid(fCalibration->GetTopPedestal(s)) || !IsValid(fCalibration->GetBotPedestal(s)) || !IsValid(fCalibration->GetTopAdcGain(s)) || !IsValid(fCalibration->GetBotAdcGain(s)) || !IsValid(fCalibration->GetTopTdcGain(s)) || !IsValid(fCalibration->GetBotTdcGain(s)) || !IsValid(fCalibration->GetTimeOffset(s))) fValidSlat[s-1] = kFALSE; } } //____________________________________________________________________ void BrTofSlewingCalModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // Per event method SetState(kEvent); Info(25,"Event","Start event, fCommitAscii = %d, fLoadAscii = %d", fCommitAscii,fLoadAscii); if (fCommitAscii || fLoadAscii) return; // ------ get the the matched tracks and hits table BrDataTable* match = inNode->GetDataTable(fMatchName.Data()); if (!match) match = outNode->GetDataTable(fMatchName.Data()); if (!match) { Info(20, "Event", "No tracks matching hits"); return; } // ------ beam beam vertex BrBbVertex* bb = (BrBbVertex*)inNode->GetObject("BB Vertex"); if (!bb) { bb = (BrBbVertex*)outNode->GetObject("BB Vertex"); } if(fUseBbVertex && !bb) { Info(10,"Event", " No beam beam stuff"); return; } // ------ select only small or big tube vertex Int_t method = 3; Double_t T0 = 0.0; Float_t Z0 = 0; if (bb){ method = bb->GetMethod(); if (fUseBbVertex && method == 3) return; T0 = bb->GetTime0(); Z0 = bb->GetZ0(); } BrDataTable* td1tab = 0; if(fUseTd1Time) td1tab = inNode->GetDataTable("FfsMatchTd1"); BrDataTable* tmrstab = 0; if(fUseTMrsTime) { tmrstab = inNode->GetDataTable("MrsFMatch"); if(Verbose()>50) inNode->ListObjects(); } Int_t ntd1=0; if (fUseTd1Time){ if(td1tab) ntd1 = td1tab->GetEntries(); else { if (DebugLevel() > 5) cout << " No TD1 hits" << endl; return; } if(Verbose()>10){ for(int k=0;k<ntd1;k++){ BrTd1TrackMatch* td1m = ( BrTd1TrackMatch*) td1tab->At(k); td1m->Print(); } } } Int_t ntmrs = 0; if (fUseTMrsTime) { if (tmrstab) ntmrs = tmrstab->GetEntries(); if (Verbose() > 10) { for(Int_t k = 0; k < ntmrs; k++) { BrTMrsFTrackMatch* tmrsm = (BrTMrsFTrackMatch*)tmrstab->At(k); tmrsm->Print(); } } } // ------ get hit table BrDataTable* hits = inNode->GetDataTable(Form("DigTof %s", GetName())); if (!hits) hits = outNode->GetDataTable(Form("DigTof %s", GetName())); if (!hits) { if (Verbose() > 20) Warning("Event", "No global tof hits at all"); return; } // ------ get track table BrDataTable* tracks = 0; if (fInFfs) tracks = inNode->GetDataTable("FfsTracks"); if (fInBfs) tracks = inNode->GetDataTable("FsTracks"); if (fInMrs) tracks = inNode->GetDataTable("MrsTracks"); if (!tracks) { if (fInFfs) tracks = outNode->GetDataTable("FfsTracks"); if (fInBfs) tracks = outNode->GetDataTable("FsTracks"); if (fInMrs) tracks = outNode->GetDataTable("MrsTracks"); } if (!tracks) { if (Verbose() > 20) Warning("Event", "No tracks at all"); return; } // ------ get Cherenkov pid table if required. BrDataTable* chkTab = 0; if(fUseChkv){ if (fInFfs) chkTab = inNode->GetDataTable("C1Pid"); if (fInBfs) chkTab = inNode->GetDataTable("RichPid"); if (!chkTab) { if (fInFfs) chkTab = outNode->GetDataTable("C1Pid"); if (fInBfs) chkTab = outNode->GetDataTable("RichPid"); } if (!chkTab) { if (Verbose() > 20) Warning("Event", "No cherenkov PID at all"); return; } } if (Verbose() > 15){ cout << " ----- Slewing correction module " << endl << " Number of tracks : " << tracks->GetEntries() << endl << " Number of tof digits : " << hits->GetEntries() << endl << " Number of matching pairs : " << match->GetEntries() << endl; if(fUseChkv) cout << " Number of Cherenkov Pid : " << chkTab->GetEntries() << endl; } // -------------------------------------------------------- // -------------------------------------------------------- Int_t nhit = hits->GetEntries(); // number of tof hits in event Int_t ntrk = tracks->GetEntries(); // number of tracks in event Int_t nchk = 0; if(fUseChkv) nchk = chkTab->GetEntries(); // number of chk pid elements Int_t nmch = match->GetEntries(); // number matching pairs vector <BrTofDig*> tofHit(nmch); // array of good hits vector <TObject*> track(nmch); // array of good tracks vector <BrPid*> chkPid(nmch); // pid object vector <BrTd1TrackMatch*> td1hit(nmch); //array of possible Td1Hits vector <BrTMrsFTrackMatch*> tmrshit(nmch); //array of possible TMrsFhits // we are now sure that the hit will correspond to a valid track // ----- loop over matches to pick up right objects for (Int_t m = 0; m < nmch; m++) { BrTofTrackMatch* pair = (BrTofTrackMatch*)match->At(m); tofHit[m] = 0; track [m] = 0; chkPid[m] = 0; // now select right hit thanks to the match object for (Int_t h = 0; h < nhit; h++) { BrTofDig* hit = (BrTofDig*)hits->At(h); if (hit->GetSlatno() != pair->GetHitId()) continue; tofHit[m] = hit; break; } // now select right track proj thanks to the match object for (Int_t t = 0; t < ntrk; t++) { TObject* trk = tracks->At(t); if (trk->GetUniqueID() != pair->GetTrackId()) continue; track[m] = trk; break; } if(fUseChkv){ // now select right pid objects for (Int_t p = 0; p < nchk; p++) { BrPid* pid = (BrPid*)chkTab->At(p); if (pid->GetTrackId() != pair->GetTrackId()) continue; chkPid[m] = pid; break; } } // Lookup the Td1Hit corresponding to this track // Note: Why is this loop backwards (copied from tof matches) if(fUseTd1Time){ if (td1tab){ for (Int_t p = ntd1-1; p >=0; p--) { BrTd1TrackMatch* chk = (BrTd1TrackMatch*)td1tab->At(p); if (chk->GetTrackId() != pair->GetTrackId()) continue; td1hit[m] = chk; break; } } } //if(fUseTd1Time) // // Lookup the TMrsFHit corresponding to this track // if(fUseTMrsTime){ if(tmrstab){ for (Int_t p = ntmrs-1; p >=0; p--) { BrTMrsFTrackMatch* chk = (BrTMrsFTrackMatch*)tmrstab->At(p); if (chk->GetTrackId() != pair->GetTrackId()) continue; tmrshit[m] = chk; } } } //if(fUseTMrsTime) } //------------------------------------------ // can now proceed with calibration //------------------------------------------ for(Int_t h = 0; h < nmch; h++) { BrTofTrackMatch* pair = (BrTofTrackMatch*)match->At(h); if (fValidSlat[pair->GetHitId() - 1] == kFALSE) continue; // ------ check if we have a valid PID in cherenkov det. if (fUseChkv && !chkPid[h]) continue; // --- track info Float_t length; Float_t moverp; Float_t momentum; Float_t chksig = 0; BrVector3D trkVtx; BrVector3D proj; // ------ if in FFS: check D1 swim status and C1 signal if (fInFfs) { BrFfsTrack* t = (BrFfsTrack*)track[h]; if (!t->GetD1SwimStatus()) continue; proj = t->GetProjOnTof(); length = t->GetPathLength(); momentum = t->GetMomentum(); moverp = 0.1395679/t->GetMomentum(); trkVtx = t->GetTrackVertex(); if (fUseChkv) { chksig = ((BrC1Pid*)chkPid[h])->GetBlobEnergy(); if (chksig < fC1Threshold) continue; } } if (fInMrs) { BrMrsTrack* t = (BrMrsTrack*)track[h]; length = t->GetPathLength(); momentum = t->GetMomentum(); moverp = 0.1395679/t->GetMomentum(); trkVtx = t->GetTrackVertex(); Int_t panel = pair->GetPanelId(); fPanelVol[panel]->GlobalToLocal(t->GetProjOnTof(), proj, 0); if(t->GetStatus() !=1) continue; } if (fInBfs) { BrFsTrack* t = (BrFsTrack*)track[h]; if (!t->GetD1SwimStatus()) continue; momentum = t->GetMomentum(); proj = t->GetProjOnTof2(); length = t->GetPathLength(); moverp = 0.1395679/t->GetMomentum(); trkVtx = t->GetTrackVertex(); if(fUseChkv){ chksig = ((BrRichPid*)chkPid[h])->GetMass(); if (TMath::Abs(chksig - fRichPionMass) > fRichMassCut) continue; } } if(fInMrs){ Int_t tmrsslat = 0; if (fUseTMrsTime) if (tmrshit[h]) { T0 = tmrshit[h]->GetTime(); length = ((BrMrsTrack*)track[h])->GetPartialPath() + tmrshit[h]->GetLength(); tmrsslat = 1; } else continue; if(TMath::Abs(momentum) > fMomentumCut) continue; // check track vertex (want only primary tracks) // (only if vertex from Bb available. // Float_t x = (trkVtx(0) - fTrkOffset[0])/fTrkCut[0]; Float_t y = (trkVtx(1) - fTrkOffset[1])/fTrkCut[1]; Float_t z = (trkVtx(2) - Z0 - fTrkOffset[2])/fTrkCut[2]; if (fUseBbVertex && (x*x + y*y + z*z > 1)) continue; } // ------ real calibration starts now Int_t slat = tofHit[h]->GetSlatno(); Float_t tped = fCalibration->GetTopPedestal(slat); Float_t bped = fCalibration->GetBotPedestal(slat); Float_t tag0 = fCalibration->GetTopAdcGain(slat); Float_t bag0 = fCalibration->GetBotAdcGain(slat); Float_t ttg = fCalibration->GetTopTdcGain(slat); Float_t btg = fCalibration->GetBotTdcGain(slat); Double_t offt =fCalibration->GetTimeOffset(slat); Double_t tadc = tofHit[h]->GetAdcUp() - tped; Double_t tener = tadc/ tag0; Double_t badc = tofHit[h]->GetAdcDown() - bped; Double_t bener = badc/bag0; if (tener < fEnergyThreshold || bener < fEnergyThreshold || tener > fMaxTofDE || bener > fMaxTofDE) continue; Double_t tTime = tofHit[h]->GetTdcUp()*ttg; Double_t bTime = tofHit[h]->GetTdcDown()*btg; // Calculate the theoretical TOF from beta // beta = p/sqrt(p^2 + m^2) with // p = momentum, m = pion mass. This means that their // tof will be aligned with t = L/(c beta), L being the path length // Either fill histograms to slewfit or fill intermediate TobjArrays // that holds all the relevant data, for fitting. // // Get the T0 and track lenght relevant for the start time detector // Double_t ffsTrackLength = 0; Double_t bfsTrackLength = 0; Double_t td1TrackLength = 0; Int_t td1slat=0; if (fUseTd1Time) if (td1hit[h]) { T0 = td1hit[h]->GetTime(); if(fInFfs) { length = ((BrFfsTrack*)track[h])->GetPartialPath() + td1hit[h]->GetLength(); td1slat = td1hit[h]->GetSlat(); } else if(fInBfs) { BrFfsTrack *ffsTrack = ((BrFsTrack*)track[h])->GetFfsTrack(); BrBfsTrack *bfsTrack = ((BrFsTrack*)track[h])->GetBfsTrack(); //We have to do following this way because there are some events //in which there is no ffsTrack in the BrFsTrack. I don't understand //that. Needs investigation. KH 25-APR-2002. In any case, //don't set td1slat if either track is missing. if(ffsTrack) ffsTrackLength = ffsTrack->GetPartialPath(); if(bfsTrack) bfsTrackLength = bfsTrack->GetPathLength(); td1TrackLength = td1hit[h]->GetLength(); length = ffsTrackLength + bfsTrackLength + td1TrackLength; if(ffsTrack && bfsTrack) td1slat = td1hit[h]->GetSlat(); } } Double_t tTheor = length/BrUnits::c_light * TMath::Sqrt(moverp*moverp+1); Float_t beta = 1.0 / TMath::Sqrt(moverp*moverp+1); Info(10,"Event","slat = %d, tTheor = %f, tMeasured = %f, offt = %f", slat,tTheor,0.5*(tTime+bTime)-T0-offt,offt); // // load eventdata that are used at end for the fits. // if(fUseMinuit){ // Check time is near predicted time. // If so increment table. // Float_t tMeasured = 0.5*(tTime+bTime)-T0-offt; if(TMath::Abs(tMeasured-tTheor)<fTimeCut) { BrTofSlewHit* obj = new BrTofSlewHit(tTime, bTime, tadc, badc, T0, length, beta); fSlewHitsTable[slat-1].Add(obj); if(HistOn()){ fTDiff->Fill(slat, tMeasured-tTheor); fBOff[slat-1]->Fill(badc, tMeasured - tTheor); fTOff[slat-1]->Fill(tadc, tMeasured - tTheor); } } } else { // Top calibration Double_t tnc = tTime - T0 - offt; // time not calibrated fTOff[slat-1]->Fill(tener, tnc - tTheor); if(TMath::Abs(tnc-tTheor)<1.0) fTSlew[slat-1]->Fill(tener, tnc - tTheor); // Bot calibration tnc = bTime - T0 -offt; fBOff[slat-1]->Fill(bener, tnc - tTheor); if(TMath::Abs(tnc-tTheor)<1.0) fBSlew[slat-1]->Fill(bener, tnc - tTheor); } } } //____________________________________________________________________ void BrTofSlewingCalModule::Finish() { //----------------------------------------------------------- // project 2d histos to profiles // fit with a + b / dE (older analyses showed that // this function gave the smallest errors on the parameters) //----------------------------------------------------------- // Job-level finalisation SetState(kFinish); // if load ascii mode if (fLoadAscii) return; // if commit mode if (fCommitAscii) { ReadAscii(); return; } TDirectory* saveDir = gDirectory; Int_t nSlats = fParamsTof->GetNoSlats(); if(fUseMinuit){ vector <Double_t> AUp(nSlats); vector <Double_t> ADown(nSlats); vector <Double_t> PowerUp(nSlats); vector <Double_t> PowerDown(nSlats); vector <Double_t> Dt(nSlats); // // Setup Minuit to be able to fit up to 5 parameters // Loop over relevant slats, discard those with bad calibration ?! // fMinuit->SetFCN(fcn2); Int_t nSlats = fParamsTof->GetNoSlats(); Double_t dfit[5]; Double_t dfitError[5]; for(Int_t slat=1; slat<=nSlats; slat++){ FitCurrentSlat(slat, dfit, dfitError); AUp[slat-1] = dfit[0]; ADown[slat-1] = dfit[2]; PowerUp[slat-1]=dfit[1]; PowerDown[slat-1]=dfit[3]; Dt[slat-1]=dfit[4]; } // Fill the calibration array with found fits // This is done quite similar to the previous method. // for(Int_t slat=0; slat<nSlats; slat++){ Double_t tpar1, tpar2, bpar1, bpar2; // check if fit failed tpar2 = AUp[slat]; bpar2 = ADown[slat]; if (tpar2==40 || bpar2==40) { tpar2 = BrTofCalibration::kCalException; bpar2 = BrTofCalibration::kCalException; tpar1 = BrTofCalibration::kCalException; bpar1 = BrTofCalibration::kCalException; } else { tpar1 = -Dt[slat]; bpar1 = -Dt[slat]; } // fill calibration parameter elements Int_t i = slat+1; fCalibration->SetTopSlewPar1(i, tpar1); fCalibration->SetTopSlewPar2(i, tpar2); fCalibration->SetBotSlewPar1(i, bpar1); fCalibration->SetBotSlewPar2(i, bpar2); fTSlew1->Fill(i, tpar1); fTSlew2->Fill(i, tpar2); fBSlew1->Fill(i, bpar1); fBSlew2->Fill(i, bpar2); dfit[0] = AUp[slat]; dfit[1] = -0.5; dfit[3]= -0.5; dfit[2] = ADown[slat]; dfit[4] = Dt[slat]; FillCorrectedHistograms(i,dfit); } } // // Histogram method with normalized Adc gain. // else{ TF1* f = new TF1("slewing", "[0] + [1]/sqrt(x)", 0.8, 2); for (Int_t i = 1; i <= fParamsTof->GetNoSlats(); i++) { if (!fValidSlat[i-1]) continue; //------------ top tubes fHistDir->cd("Top"); Double_t tpar1, tpar2, bpar1, bpar2; TProfile* p = fTSlew[i-1]; f->SetParameters(5, 0.5); p->Fit("slewing", "Q0", "", fEnergyThreshold, fMaxTofDE); tpar1 = f->GetParameter(0); tpar2 = f->GetParameter(1); if (DebugLevel() > 0) cout << " top " << i << " a = " << f->GetParameter(0) << " b = " << f->GetParameter(1); //------------ bot tubes fHistDir->cd("Bot"); p = fBSlew[i-1]; f->SetParameters(5, 0.5); p->Fit("slewing", "Q0", "", fEnergyThreshold, fMaxTofDE); bpar1 = f->GetParameter(0); bpar2 = f->GetParameter(1); if (DebugLevel() > 0) cout << " --- bot " << i << " a = " << f->GetParameter(0) << " b = " << f->GetParameter(1) << endl; // ----------------------------------------------- // check if fit failed if (tpar1==5 || tpar2==0.5) { tpar1 = BrTofCalibration::kCalException; tpar2 = BrTofCalibration::kCalException; } if (bpar1==5 || bpar2==0.5) { bpar1 = BrTofCalibration::kCalException; bpar2 = BrTofCalibration::kCalException; } // fill calibration parameter elements fCalibration->SetTopSlewPar1(i, tpar1); fCalibration->SetTopSlewPar2(i, tpar2); fCalibration->SetBotSlewPar1(i, bpar1); fCalibration->SetBotSlewPar2(i, bpar2); fTSlew1->Fill(i, tpar1); fTSlew2->Fill(i, tpar2); fBSlew1->Fill(i, bpar1); fBSlew2->Fill(i, bpar2); if (DebugLevel()) cout << " Top: par1 -> " << fCalibration->GetTopSlewPar1(i) << " par2 -> " << fCalibration->GetTopSlewPar2(i) << endl << " Bot: par1 -> " << fCalibration->GetBotSlewPar1(i) << " par2 -> " << fCalibration->GetBotSlewPar2(i) << endl; } } gDirectory = saveDir; // save calibration to ascii file if (fSaveAscii) SaveAscii(); } //___________________________________________________________________________ void BrTofSlewingCalModule::FillCorrectedHistograms(Int_t slat, Double_t* par){ // // Take Hit data and apply correction to data and store // in histograms. Note the // fCurrentHitTable = &fSlewHitsTable[slat-1]; fCurrentSlat = slat; int np =fCurrentHitTable->GetEntries(); if(DebugLevel()>5){ cout << "Filling Corrected Histograms " << setw(4) << slat << setw(8) << np << endl; cout << "Params " << par[0] << setw(8) << par[2] << setw(8) << par[4] << endl; } for(Int_t i; i< np; i++) { BrTofSlewHit* hit = (BrTofSlewHit*) fCurrentHitTable->At(i); Double_t value = 0.5*(hit->fTofup + hit->fTofdown) -0.5*( par[0]/TMath::Sqrt(hit->fAdcup) + par[2]/TMath::Sqrt(hit->fAdcdown)) -fCalibration->GetTimeOffset(fCurrentSlat) -hit->fT0 + par[4]; // fitmode 1 that subtracts nominal dt change. if(fFitMode==1) value -= 0.5*(par[0]/TMath::Sqrt(fCalibration->GetBotAdcGain(fCurrentSlat)) +par[2]/TMath::Sqrt(fCalibration->GetTopAdcGain(fCurrentSlat))); Double_t tTheor = hit->fLength / hit->fBetaCalc/BrUnits::c_light; fTOffCorr[slat-1]->Fill(hit->fAdcup, value - tTheor); fBOffCorr[slat-1]->Fill(hit->fAdcdown, value - tTheor); fTDiffCorr->Fill(slat, value - tTheor); } } //____________________________________________________________________ void BrTofSlewingCalModule::FitCurrentSlat(Int_t slat, Double_t* dfit, Double_t* dfitError) { // // // set the global instance to this module, so the fitting is done with the // correct object members; // fgInstance = this; fCurrentHitTable = &fSlewHitsTable[slat-1]; fCurrentSlat = slat; if(Verbose()>2){ cout << "Slat # " << slat << "n" << " Np " << fCurrentHitTable->GetEntries() << endl; } Double_t vstart[] = { 40., -0.5, 40., -0.5, 0.0}; Double_t step[] = { 2.5, 0.05, 2.5, 0.05, 0.1}; for(int k=0;k<5;k++){ dfit[k]=vstart[k]; dfitError[k]=-1.0; } Int_t np = fCurrentHitTable->GetEntries(); if(np < 20) { cout<<"Warning: for slat "<<slat<<" Number of entries < 50; it is "; cout<<np<<endl; if(np < 50) return; } Int_t ierflg; // Equivalent Miniut command- fMinuit->mnexcm("SET ERR", arglist, 1, ierflg); // Since the virtual base clas is call cannot do specific calls !! // The virtual class is use so not to include minuit in libraries unless needed. // Double_t arglist[2]; arglist[0] = 1; fMinuit->ExecuteCommand("SET ERR", arglist ,1); arglist[0]=2000.; arglist[1]=0.1; Int_t errorFlag; // return from minuit // set initial values // They correspond approximately to the value of 0.5/sqrt(Adc/gain) // errorFlag = fMinuit->SetParameter(0, "Const(up) " , vstart[0], step[0], 4, 100); errorFlag = fMinuit->SetParameter(2, "Const(down)" , vstart[2], step[2], 4, 100); errorFlag = fMinuit->SetParameter(1, "Power(up) " , vstart[1], step[1], 0, 0); errorFlag = fMinuit->SetParameter(3, "Power(down)" , vstart[3], step[3], 0, 0); errorFlag = fMinuit->SetParameter(4, "delta(t)" , vstart[4], step[4], 0, 0); // assume a fixed power coefficient for now. // fMinuit->FixParameter(1); fMinuit->FixParameter(3); if(fVaryDt) fMinuit->ReleaseParameter(4); else fMinuit->FixParameter(4); // perform minimization // fMinuit->mnexcm("MIGRAD", arglist, 2, errorFlag); errorFlag = fMinuit->ExecuteCommand("MIGRAD", arglist ,2); fMinuit->FixParameter(2); fMinuit->ReleaseParameter(2); if(Verbose()>10){ /* arglist[0]=0; errorFlag = fMinuit->ExecuteCommand("SCAN", arglist ,1); arglist[0]=2; errorFlag = fMinuit->ExecuteCommand("SCAN", arglist ,1); */ arglist[0]=0; arglist[1]=2; errorFlag = fMinuit->ExecuteCommand("CONTOUR", arglist ,2); } // Extract and store result // How to check for successfull minimization. Char_t name[32]; Double_t low, high; fMinuit->GetParameter(0,name, dfit[0], dfitError[0], low, high); fMinuit->GetParameter(2,name, dfit[2], dfitError[2], low, high); if(fVaryDt) fMinuit->GetParameter(4,name, dfit[4], dfitError[4], low, high); else { // Determine dt from mean deviation of fitting function from expected function // Double_t mean, meanwt, chisq, chisqpdf; fcnout(dfit, mean, meanwt, chisq, chisqpdf); dfit[4] = (fNominalLength/BrUnits::c_light - mean)/meanwt; if(Verbose()>1){ cout << " Dt determination" << setw(12) <<mean << setw(12) << meanwt << setw(12) << chisq << setw(12) <<chisqpdf << endl; } } if(Verbose()>1){ cout << " Fit result " << setw(8) << dfit[0] << " +-" << dfitError[0] << endl; cout << " " << setw(8) << dfit[2] << " +-" << dfitError[2] << endl; cout << " " << setw(8) << dfit[4] << " +-" << dfitError[4] << endl; } } //____________________________________________________________________ void BrTofSlewingCalModule::fcn2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) { // static function so it can only operate on static objects, meaning // that all objects/variables must be known. Trick is to make a // static objectpointer // calculate chisquare f = BrTofSlewingCalModule::Instance()->ChiSquare(par); } //_________________________________________________________________________ Double_t BrTofSlewingCalModule::ChiSquare(Double_t *par) { // This function calculates the chisquare for all fit points // Most of these parameters are needed for Minuit- they are not for specific // use by the program. // Double_t mean = 0.0; Double_t chisq = 0.0; Double_t delta; Int_t np = fCurrentHitTable->GetEntries(); for(Int_t i=0; i<np; i++) mean +=func2(i, par); mean= mean / Double_t (np); for(Int_t i=0; i<np; i++){ Double_t delta = (mean - func2(i,par))/ fTimeResolution; chisq += delta*delta; } if(DebugLevel()>5){ cout << "param " << setw(8) << par[0] << setw(8) << par[2] << setw(8) << par[4] << endl; cout << "mean " << mean << " np=" << np << endl; cout << "chisq " << setw(8) << chisq << endl; } return chisq; } //____________________________________________________________________ void BrTofSlewingCalModule::fcnout(Double_t *par, Double_t &mean, Double_t &meanwt, Double_t &chisq, Double_t &chisqpdf){ // // Calculate final mean time, acerave tweight to get the // Assumes the the currenthitable already point to the correct slat. // Should be ok this fct is called from FitCurrentSlat. // Double_t delta; mean = 0.0; meanwt = 0.0; Int_t np = fCurrentHitTable->GetEntries(); for(Int_t i=0; i<np; i++){ BrTofSlewHit* hit = (BrTofSlewHit*) fCurrentHitTable->At(i); mean +=func2(i, par); meanwt += (hit->fBetaCalc*fNominalLength / hit->fLength); } mean = mean / Double_t (np); meanwt = meanwt / Double_t (np); chisq = 0.0; for(Int_t i=0; i<np; i++){ Double_t delta = (mean - func2(i,par))/ fTimeResolution; chisq += delta*delta; } chisqpdf = chisq/Double_t(np-2); } //____________________________________________________________________ Double_t BrTofSlewingCalModule::func2(Int_t i, Double_t *par) { // Evaluate beta(calc)/beta(measured) - <beta(calc)/beta(measured)> // for datapoint i. // BrTofSlewHit* hit = (BrTofSlewHit*) fCurrentHitTable->At(i); Double_t value = 0.5*(hit->fTofup + hit->fTofdown) -0.5*( par[0]/TMath::Sqrt(hit->fAdcup)+ par[2]/TMath::Sqrt(hit->fAdcdown)) -fCalibration->GetTimeOffset(fCurrentSlat) -hit->fT0 + par[4]; if(fFitMode==1) value -= 0.5*(par[0]/TMath::Sqrt(fCalibration->GetBotAdcGain(fCurrentSlat)) +par[2]/TMath::Sqrt(fCalibration->GetTopAdcGain(fCurrentSlat))); if(DebugLevel()>10){ cout << " hit" << setw(8)<<hit->fTofup << setw(8) << hit->fTofdown << setw(8) << hit->fAdcup<< setw(8) << hit->fAdcdown<< setw(8) << hit->fT0 <<setw(8) << hit->fBetaCalc << endl; cout << " " << setw(8)<<0.5*(hit->fTofup+hit->fTofdown)-fCalibration->GetTimeOffset(fCurrentSlat) -hit->fT0 << endl; cout << " " << setw(8) << value; } value *= (hit->fBetaCalc*fNominalLength/hit->fLength); if(DebugLevel()>10){ cout << " " << setw(8) << value << endl; } return value; } //____________________________________________________________________ void BrTofSlewingCalModule::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 << "* Slewing parameters " << 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 par1 | top par2 | bot par1 | bot par2" << endl; file << "* -------------------------------------------------------" << endl << endl; for (Int_t i = 0; i < fParamsTof->GetNoSlats(); i++) { Int_t slat = i + 1; file << setw(4) << slat << setw(15) << fCalibration->GetTopSlewPar1(slat) << setw(15) << fCalibration->GetTopSlewPar2(slat) << setw(15) << fCalibration->GetBotSlewPar1(slat) << setw(15) << fCalibration->GetBotSlewPar2(slat) << endl; } file << "* ------------------------------------" << endl << endl; } //____________________________________________________________________ void BrTofSlewingCalModule::ReadAscii() { BrTofCalModule::ReadAscii(); ifstream file(fCalibFile.Data(), ios::in); if (!file) { Stop("ReadAscii", "File %s was not found", fCalibFile.Data()); return; } Float_t tp1, tp2, bp1, bp2; 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 >> tp1 >> tp2 >> bp1 >> bp2; fCalibration->SetTopSlewPar1(slat, tp1); fCalibration->SetTopSlewPar2(slat, tp2); fCalibration->SetBotSlewPar1(slat, bp1); fCalibration->SetBotSlewPar2(slat, bp2); if (DebugLevel() > 5) cout << setw(4) << slat << setw(12) << tp1 << setw(12) << tp2 << setw(12) << bp1 << setw(12) << bp2 << endl; } fCalibration->SetComment("topSlewPar1", fComment.Data()); fCalibration->SetComment("topSlewPar2", fComment.Data()); fCalibration->SetComment("botSlewPar1", fComment.Data()); fCalibration->SetComment("botSlewPar2", fComment.Data()); } //____________________________________________________________________ void BrTofSlewingCalModule::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: hagel $" << endl << " $Date: 2002/08/06 19:04:30 $" << endl << " $Revision: 1.16 $ " << endl << endl << "-------------------------------------------------" << endl; } //____________________________________________________________________ // // $Log: BrTofSlewingCalModule.cxx,v $ // Revision 1.16 2002/08/06 19:04:30 hagel // Implement slewing correction calibration for TOFW // // Revision 1.15 2002/07/01 16:01:38 hagel // Implement td1 timing ala pp for H2 calibration // // Revision 1.14 2002/06/19 17:34:27 videbaek // Fix momentum cut (make it an absolute value - MRS // increase no bins in timediff histogarms to 40 psec/bin // // Revision 1.13 2002/06/14 17:48:08 videbaek // Fixed a couple of instanmces of wrongly used no slats (in new minuit method). // // Revision 1.12 2002/06/13 19:24:15 videbaek // mny modification to timeoffset and slewing modules. // a) Offset module modified to deal with pp running and Td1, TMrsF start counters. // b) The slewing and offset mdoule now has similar cuts in pid, vertex cuts // and momentmu cut. // c) The mrs code had a momentum cut introduced to not contaminate dt spectra with // the kaon's and protons. // d) The slewing modules have added the method to fit the hits using minuit rather than // 2D histograms. the old method is kept for consistency (at least a while). As for // all slwing data correction rather large statistics is needed. // e) Note that when using slewing the timeoffsets are still important!!! the first term in the // slweing is only a relative term taking into account the offsets created by the 1/sqrt(A) terms. // This implies that the same offsets are valid w and /wo slewing being done. // // Revision 1.11 2002/04/15 17:15:56 ouerdane // added SetMaxEnergy(max tof de) // // Revision 1.10 2002/04/09 02:03:10 ouerdane // use BrFsTrack instead of BrBfsTrack for H2 calibration // // Revision 1.9 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.8 2002/03/20 19:37:25 videbaek // Added possibility to calibrate for pp using the TD1 counters rather than // BB. Also falg possibility to require or not Cherenkovs. This may be useful // for calibrations in MRS. // // Revision 1.7 2001/12/20 15:30:44 ouerdane // Modified slewing correction, needs now the Rich or C1 pid objects to select only pions (cf calib/tof/TofSlewing.C // // Revision 1.6 2001/11/05 06:59:44 ouerdane // changed to deal with new track classes and fixed some bugs // // Revision 1.5 2001/10/20 00:28:16 ouerdane // increased verbosity level to > 10 // // Revision 1.4 2001/10/08 11:28:10 cholm // Changed to use new DB access classes // // Revision 1.3 2001/10/02 01:53:28 ouerdane // Added SetSaveAscii, SetLoadAscii, SetCommitAscii and updated the way parameters are tagged for Use // // Revision 1.2 2001/07/31 09:22:00 ouerdane // Removed all references to a TList member, replaced // by histogram members, declare ntuple if fNtuple is true (set // via SetNtuple) // // Revision 1.1 2001/06/29 13:52:25 cholm // Imported TOF classes from Djamel // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|