![]() |
![]() |
#include <CalibMuon/DTCalibration/plugins/DTT0CalibrationNew.h>
Public Member Functions | |
void | analyze (const edm::Event &event, const edm::EventSetup &eventSetup) |
Fill the maps with t0 (by channel). | |
DTT0CalibrationNew (const edm::ParameterSet &pset) | |
Constructor. | |
void | endJob () |
Compute the mean and the RMS of the t0 from the maps and write them to the DB with channel granularity. | |
virtual | ~DTT0CalibrationNew () |
Destructor. | |
Private Member Functions | |
std::string | getHistoName (const DTLayerId &lId) const |
std::string | getHistoName (const DTWireId &wId) const |
Private Attributes | |
std::vector< std::string > | cellsWithHistos |
bool | debug |
std::string | digiLabel |
edm::ESHandle< DTGeometry > | dtGeom |
unsigned int | eventsForLayerT0 |
unsigned int | eventsForWireT0 |
TH1D * | hT0SectorHisto |
std::map< DTWireId, double > | mK |
std::map< DTWireId, double > | mK_ref |
std::map< DTWireId, int > | nDigiPerWire |
std::map< DTWireId, int > | nDigiPerWire_ref |
unsigned int | nevents |
std::map< DTWireId, double > | qK |
unsigned int | rejectDigiFromPeak |
unsigned int | retryForLayerT0 |
int | selSector |
int | selWheel |
TSpectrum * | spectrum |
std::map< DTWireId, double > | theAbsoluteT0PerWire |
std::string | theCalibSector |
std::string | theCalibWheel |
std::map< DTChamberId, int > | theCountT0ByChamber |
TFile * | theFile |
std::map< DTLayerId, TH1I * > | theHistoLayerMap |
std::map< DTWireId, TH1I * > | theHistoWireMap |
TFile * | theOutputFile |
std::map< DTChamberId, double > | theRefT0ByChamber |
std::map< DTWireId, double > | theRelativeT0PerWire |
std::map< std::string, double > | theSigmaT0LayerMap |
std::map< DTWireId, double > | theSigmaT0PerWire |
std::map< DTChamberId, double > | theSumT0ByChamber |
std::map< std::string, double > | theT0LayerMap |
std::map< DTLayerId, double > | theTPPeakMap |
unsigned int | timeBoxWidth |
double | tpPeakWidth |
double | tpPeakWidthPerLayer |
std::vector< DTWireId > | wireIdWithHistos |
Those values are written in the DB with cell granularity. The mean value for each channel is normalized to a reference time common to all the sector. The t0 of wires in odd layers are corrected for the relative difference between odd and even layers
Definition at line 31 of file DTT0CalibrationNew.h.
DTT0CalibrationNew::DTT0CalibrationNew | ( | const edm::ParameterSet & | pset | ) |
Constructor.
Definition at line 36 of file DTT0CalibrationNew.cc.
References cellsWithHistos, GenMuonPlsPt100GeV_cfg::cout, debug, digiLabel, lat::endl(), eventsForLayerT0, eventsForWireT0, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), hT0SectorHisto, nevents, rejectDigiFromPeak, retryForLayerT0, selSector, selWheel, sl, spectrum, theCalibSector, theCalibWheel, theFile, timeBoxWidth, tpPeakWidth, tpPeakWidthPerLayer, muonGeometry::wheel, and wireIdWithHistos.
00036 { 00037 // Get the debug parameter for verbose output 00038 debug = pset.getUntrackedParameter<bool>("debug"); 00039 if(debug) 00040 cout << "[DTT0CalibrationNew]Constructor called!" << endl; 00041 00042 // Get the label to retrieve digis from the event 00043 digiLabel = pset.getUntrackedParameter<string>("digiLabel"); 00044 00045 // The root file which contain the histos per layer 00046 string rootFileName = pset.getUntrackedParameter<string>("rootFileName","DTT0PerLayer.root"); 00047 theFile = new TFile(rootFileName.c_str(), "RECREATE"); 00048 00049 theCalibWheel = pset.getUntrackedParameter<string>("calibWheel", "All"); //FIXME amke a vector of integer instead of a string 00050 if(theCalibWheel != "All") { 00051 stringstream linestr; 00052 int selWheel; 00053 linestr << theCalibWheel; 00054 linestr >> selWheel; 00055 cout << "[DTT0CalibrationNewPerLayer] chosen wheel " << selWheel << endl; 00056 } 00057 00058 // Sector/s to calibrate 00059 theCalibSector = pset.getUntrackedParameter<string>("calibSector", "All"); //FIXME amke a vector of integer instead of a string 00060 if(theCalibSector != "All") { 00061 stringstream linestr; 00062 int selSector; 00063 linestr << theCalibSector; 00064 linestr >> selSector; 00065 cout << "[DTT0CalibrationNewPerLayer] chosen sector " << selSector << endl; 00066 } 00067 00068 vector<string> defaultCell; 00069 defaultCell.push_back("None"); 00070 cellsWithHistos = pset.getUntrackedParameter<vector<string> >("cellsWithHisto", defaultCell); 00071 for(vector<string>::const_iterator cell = cellsWithHistos.begin(); cell != cellsWithHistos.end(); cell++){ 00072 if((*cell) != "None"){ 00073 stringstream linestr; 00074 int wheel,sector,station,sl,layer,wire; 00075 linestr << (*cell); 00076 linestr >> wheel >> sector >> station >> sl >> layer >> wire; 00077 wireIdWithHistos.push_back(DTWireId(wheel,station,sector,sl,layer,wire)); 00078 } 00079 } 00080 00081 hT0SectorHisto=0; 00082 00083 nevents=0; 00084 eventsForLayerT0 = pset.getParameter<unsigned int>("eventsForLayerT0"); 00085 eventsForWireT0 = pset.getParameter<unsigned int>("eventsForWireT0"); 00086 tpPeakWidth = pset.getParameter<double>("tpPeakWidth"); 00087 tpPeakWidthPerLayer = pset.getParameter<double>("tpPeakWidthPerLayer"); 00088 timeBoxWidth = pset.getParameter<unsigned int>("timeBoxWidth"); 00089 rejectDigiFromPeak = pset.getParameter<unsigned int>("rejectDigiFromPeak"); 00090 00091 spectrum = new TSpectrum(5); 00092 retryForLayerT0 = 0; 00093 }
DTT0CalibrationNew::~DTT0CalibrationNew | ( | ) | [virtual] |
Destructor.
Definition at line 96 of file DTT0CalibrationNew.cc.
References GenMuonPlsPt100GeV_cfg::cout, debug, lat::endl(), spectrum, and theFile.
00096 { 00097 if(debug) 00098 cout << "[DTT0CalibrationNew]Destructor called!" << endl; 00099 00100 delete spectrum; 00101 theFile->Close(); 00102 }
void DTT0CalibrationNew::analyze | ( | const edm::Event & | event, | |
const edm::EventSetup & | eventSetup | |||
) | [virtual] |
Fill the maps with t0 (by channel).
Perform the real analysis.
Implements edm::EDAnalyzer.
Definition at line 105 of file DTT0CalibrationNew.cc.
References funct::abs(), DTSuperLayerId::chamberId(), GenMuonPlsPt100GeV_cfg::cout, debug, digiLabel, dtGeom, lat::endl(), edm::EventID::event(), eventsForLayerT0, eventsForWireT0, find(), edm::EventSetup::get(), getHistoName(), edm::Event::id(), it, mean(), mK, mK_ref, nDigiPerWire, nDigiPerWire_ref, nevents, qK, rejectDigiFromPeak, edm::EventID::run(), DTChamberId::sector(), selSector, selWheel, python::multivaluedict::sort(), spectrum, DTLayerId::superlayerId(), theAbsoluteT0PerWire, theCalibSector, theCalibWheel, theCountT0ByChamber, theFile, theHistoLayerMap, theHistoWireMap, theSumT0ByChamber, theTPPeakMap, timeBoxWidth, tpPeakWidth, tpPeakWidthPerLayer, DTChamberId::wheel(), and wireIdWithHistos.
00105 { 00106 if(debug || event.id().event() % 500==0) 00107 cout << "--- [DTT0CalibrationNew] Analysing Event: #Run: " << event.id().run() 00108 << " #Event: " << event.id().event() << endl; 00109 nevents++; 00110 00111 // Get the digis from the event 00112 Handle<DTDigiCollection> digis; 00113 event.getByLabel(digiLabel, digis); 00114 00115 // Get the DT Geometry 00116 eventSetup.get<MuonGeometryRecord>().get(dtGeom); 00117 00118 // Get ttrig DB 00119 edm::ESHandle<DTTtrig> tTrigMap; 00120 eventSetup.get<DTTtrigRcd>().get(tTrigMap); 00121 00122 // Iterate through all digi collections ordered by LayerId 00123 DTDigiCollection::DigiRangeIterator dtLayerIt; 00124 for (dtLayerIt = digis->begin(); 00125 dtLayerIt != digis->end(); 00126 ++dtLayerIt){ 00127 // Get the iterators over the digis associated with this LayerId 00128 const DTDigiCollection::Range& digiRange = (*dtLayerIt).second; 00129 00130 // Get the layerId 00131 const DTLayerId layerId = (*dtLayerIt).first; //FIXME: check to be in the right sector 00132 const DTChamberId chamberId = layerId.superlayerId().chamberId(); 00133 00134 if((theCalibWheel != "All") && (layerId.superlayerId().chamberId().wheel() != selWheel)) 00135 continue; 00136 if((theCalibSector != "All") && (layerId.superlayerId().chamberId().sector() != selSector)) 00137 continue; 00138 00139 //if(debug) { 00140 // cout << "Layer " << layerId<<" with "<<distance(digiRange.first, digiRange.second)<<" digi"<<endl; 00141 //} 00142 00143 float tTrig,tTrigRMS; 00144 tTrigMap->slTtrig(layerId.superlayerId(), tTrig, tTrigRMS); 00145 if(debug&&(nevents <= 1)){ 00146 cout << " Superlayer: " << layerId.superlayerId() << endl 00147 << " tTrig,tTrigRMS= " << tTrig << ", " << tTrigRMS << endl; 00148 } 00149 00150 // Loop over all digis in the given layer 00151 for (DTDigiCollection::const_iterator digi = digiRange.first; 00152 digi != digiRange.second; 00153 digi++) { 00154 double t0 = (*digi).countsTDC(); 00155 00156 //Use first bunch of events to fill t0 per layer 00157 if(nevents < eventsForLayerT0){ 00158 //Get the per-layer histo from the map 00159 TH1I *hT0LayerHisto = theHistoLayerMap[layerId]; 00160 //If it doesn't exist, book it 00161 if(hT0LayerHisto == 0){ 00162 theFile->cd(); 00163 float hT0Min = tTrig - 2*tTrigRMS; 00164 float hT0Max = hT0Min + timeBoxWidth; 00165 hT0LayerHisto = new TH1I(getHistoName(layerId).c_str(), 00166 "T0 from pulses by layer (TDC counts, 1 TDC count = 0.781 ns)", 00167 timeBoxWidth,hT0Min,hT0Max); 00168 if(debug) 00169 cout << " New T0 per Layer Histo: " << hT0LayerHisto->GetName() << endl; 00170 theHistoLayerMap[layerId] = hT0LayerHisto; 00171 } 00172 00173 //Fill the histos 00174 theFile->cd(); 00175 if(hT0LayerHisto != 0) { 00176 // if(debug) 00177 // cout<<"Filling histo "<<hT0LayerHisto->GetName()<<" with digi "<<t0<<" TDC counts"<<endl; 00178 hT0LayerHisto->Fill(t0); 00179 } 00180 } 00181 00182 //Use all the remaining events to compute t0 per wire 00183 if(nevents>eventsForLayerT0){ 00184 // Get the wireId 00185 const DTWireId wireId(layerId, (*digi).wire()); 00186 if(debug) { 00187 cout << " Wire: " << wireId << endl 00188 << " time (TDC counts): " << (*digi).countsTDC()<< endl; 00189 } 00190 00191 //Fill the histos per wire for the chosen cells 00192 vector<DTWireId>::iterator it = find(wireIdWithHistos.begin(),wireIdWithHistos.end(),wireId); 00193 if (it!=wireIdWithHistos.end()){ 00194 //Get the per-wire histo from the map 00195 TH1I *hT0WireHisto = theHistoWireMap[wireId]; 00196 //If it doesn't exist, book it 00197 if(hT0WireHisto == 0){ 00198 theFile->cd(); 00199 hT0WireHisto = new TH1I(getHistoName(wireId).c_str(),"T0 from pulses by wire (TDC counts, 1 TDC count = 0.781 ns)",7000,0,7000); 00200 //hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())-100, 00201 //hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())+100); 00202 if(debug) 00203 cout << " New T0 per wire Histo: " << hT0WireHisto->GetName() << endl; 00204 theHistoWireMap[wireId] = hT0WireHisto; 00205 } 00206 //Fill the histos 00207 theFile->cd(); 00208 if(hT0WireHisto != 0) { 00209 //if(debug) 00210 // cout<<"Filling histo "<<hT0WireHisto->GetName()<<" with digi "<<t0<<" TDC counts"<<endl; 00211 hT0WireHisto->Fill(t0); 00212 } 00213 } 00214 00215 //Check the tzero has reasonable value 00216 //float hT0Min = tTrig - 2*tTrigRMS; 00217 //float hT0Max = hT0Min + timeBoxWidth; 00218 /*if(abs(hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin()) - t0) > rejectDigiFromPeak){ 00219 if(debug) 00220 cout<<"digi skipped because t0 per sector "<<hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())<<endl; 00221 continue; 00222 }*/ 00223 /*if((t0 < hT0Min)||(t0 > hT0Max)){ 00224 if(debug) 00225 cout<<"digi skipped because t0 outside of interval (" << hT0Min << "," << hT0Max << ")" <<endl; 00226 continue; 00227 }*/ 00228 //Select per layer 00229 if(fabs(theTPPeakMap[layerId] - t0) > rejectDigiFromPeak){ 00230 if(debug) 00231 cout<<"digi skipped because t0 too far from peak " << theTPPeakMap[layerId] << endl; 00232 continue; 00233 } 00234 00235 //Find to ref. per chamber 00236 theSumT0ByChamber[chamberId] = theSumT0ByChamber[chamberId] + t0; 00237 theCountT0ByChamber[chamberId]++; 00238 00239 //Use second bunch of events to compute a t0 reference per wire 00240 if(nevents< (eventsForLayerT0 + eventsForWireT0)){ 00241 if(!nDigiPerWire_ref[wireId]){ 00242 mK_ref[wireId] = 0; 00243 } 00244 nDigiPerWire_ref[wireId] = nDigiPerWire_ref[wireId] + 1; 00245 mK_ref[wireId] = mK_ref[wireId] + (t0-mK_ref[wireId])/nDigiPerWire_ref[wireId]; 00246 } 00247 //Use last all the remaining events to compute the mean and sigma t0 per wire 00248 else if(nevents>(eventsForLayerT0 + eventsForWireT0)){ 00249 if(abs(t0-mK_ref[wireId]) > tpPeakWidth) 00250 continue; 00251 if(!nDigiPerWire[wireId]){ 00252 theAbsoluteT0PerWire[wireId] = 0; 00253 qK[wireId] = 0; 00254 mK[wireId] = 0; 00255 } 00256 nDigiPerWire[wireId] = nDigiPerWire[wireId] + 1; 00257 theAbsoluteT0PerWire[wireId] = theAbsoluteT0PerWire[wireId] + t0; 00258 //theSigmaT0PerWire[wireId] = theSigmaT0PerWire[wireId] + (t0*t0); 00259 qK[wireId] = qK[wireId] + ((nDigiPerWire[wireId]-1)*(t0-mK[wireId])*(t0-mK[wireId])/nDigiPerWire[wireId]); 00260 mK[wireId] = mK[wireId] + (t0-mK[wireId])/nDigiPerWire[wireId]; 00261 } 00262 }//end if(nevents>1000) 00263 }//end loop on digi 00264 }//end loop on layer 00265 00266 //Use the t0 per layer histos to have an indication about the t0 position 00267 if(nevents == eventsForLayerT0){ 00268 bool increaseEvtsForLayerT0 = false; 00269 for(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin(); 00270 lHisto != theHistoLayerMap.end(); 00271 lHisto++){ 00272 if(debug) 00273 cout<<"Reading histogram "<<(*lHisto).second->GetName()<<" with mean "<<(*lHisto).second->GetMean()<<" and RMS "<<(*lHisto).second->GetRMS() << endl; 00274 00275 //Find peaks 00276 //int npeaks = spectrum->Search((*lHisto).second,0.5,"goff"); 00277 //int npeaks = spectrum->Search((*lHisto).second,(tpPeakWidthPerLayer/2.),"goff",0.3); 00278 int npeaks = spectrum->Search((*lHisto).second,(tpPeakWidthPerLayer/2.),"",0.3); 00279 00280 float *peaks = spectrum->GetPositionX(); 00281 //Put in a std::vector<float> 00282 vector<float> peakMeans(peaks,peaks + npeaks); 00283 //Sort the peaks in ascending order 00284 sort(peakMeans.begin(),peakMeans.end()); 00285 00286 //Find best peak -- preliminary criteria: find peak closest to center of time box 00287 float tTrig,tTrigRMS; 00288 tTrigMap->slTtrig((*lHisto).first.superlayerId(), tTrig, tTrigRMS); 00289 float timeBoxCenter = (2*tTrig + (float)timeBoxWidth)/2.; 00290 float hMin = (*lHisto).second->GetXaxis()->GetXmin(); 00291 float hMax = (*lHisto).second->GetXaxis()->GetXmax(); 00292 vector<float>::const_iterator tpPeak = peakMeans.end(); 00293 for(vector<float>::const_iterator it = peakMeans.begin(); it != peakMeans.end(); ++it){ 00294 float mean = *it; 00295 00296 int bin = (*lHisto).second->GetXaxis()->FindBin(mean); 00297 float yp = (*lHisto).second->GetBinContent(bin); 00298 if(debug) cout << "Peak : (" << mean << "," << yp << ")" << endl; 00299 00300 //Find RMS 00301 float previous_peak = (it == peakMeans.begin())?hMin:*(it - 1); 00302 float next_peak = (it == (peakMeans.end()-1))?hMax:*(it + 1); 00303 00304 float rangemin = mean - (mean - previous_peak)/8.; 00305 float rangemax = mean + (next_peak - mean)/8.; 00306 int binmin = (*lHisto).second->GetXaxis()->FindBin(rangemin); 00307 int binmax = (*lHisto).second->GetXaxis()->FindBin(rangemax); 00308 (*lHisto).second->GetXaxis()->SetRange(binmin,binmax); 00309 //RMS estimate 00310 float rms_seed = (*lHisto).second->GetRMS(); 00311 00312 /*rangemin = mean - 2*rms_seed; 00313 rangemax = mean + 2*rms_seed; 00314 if(debug) cout << "Seed for RMS, Fit min, Fit max: " << rms_seed << ", " << rangemin << ", " << rangemax << endl; 00315 //Fit to gaussian 00316 string funcname("fitFcn_"); 00317 funcname += (*lHisto).second->GetName(); 00318 if(debug) cout << "Fitting function " << funcname << endl; 00319 TF1* func = new TF1(funcname.c_str(),"gaus",rangemin,rangemax); 00320 func->SetParameters(yp,mean,rms_seed); 00321 (*lHisto).second->Fit(func,"Q","",rangemin,rangemax); 00322 00323 float fitconst = func->GetParameter(0); 00324 float fitmean = func->GetParameter(1); 00325 float fitrms = func->GetParameter(2); 00326 float chisquare = func->GetChisquare()/func->GetNDF(); 00327 if(debug) cout << "Gaussian fit constant,mean,RMS,chi2= " << fitconst << ", " << fitmean << ", " << fitrms << ", " << chisquare << endl;*/ 00328 00329 //Reject peaks with RMS larger than specified 00330 //if(fitrms > tpPeakWidth) continue; 00331 if(rms_seed > tpPeakWidthPerLayer) continue; 00332 00333 if(fabs(mean - timeBoxCenter) < fabs(*tpPeak - timeBoxCenter)) tpPeak = it; 00334 } 00335 //Didn't find peak 00336 /*if(tpPeak == peakMeans.end()){ 00337 if(retryForLayerT0 < 2){ 00338 increaseEvtsForLayerT0 = true; 00339 retryForLayerT0++; 00340 break; 00341 } 00342 }*/ 00343 00344 float selPeak = (tpPeak != peakMeans.end())?*tpPeak:(*lHisto).second->GetBinCenter((*lHisto).second->GetMaximumBin()); 00345 if(debug) cout << "Peak selected at " << selPeak << endl; 00346 00347 theTPPeakMap[(*lHisto).first] = selPeak; 00348 00349 //Take the mean as a first t0 estimation 00350 /*if((*lHisto).second->GetRMS() < tpPeakWidth){ 00351 if(hT0SectorHisto == 0){ 00352 hT0SectorHisto = new TH1D("hT0AllLayerOfSector","T0 from pulses per layer in sector", 00353 //20, (*lHisto).second->GetMean()-100, (*lHisto).second->GetMean()+100); 00354 700, 0, 7000); 00355 //300,3300,3600); 00356 } 00357 if(debug) 00358 cout<<" accepted"<<endl; 00359 //TH1I* aux_histo = (*lHisto).second; 00360 //aux_histo->GetXaxis()->SetRangeUser(3300,3600); 00361 hT0SectorHisto->Fill((*lHisto).second->GetMean()); 00362 //hT0SectorHisto->Fill(aux_histo->GetMean()); 00363 } 00364 //Take the mean of noise + 400ns as a first t0 estimation 00365 //if((*lHisto).second->GetRMS()>10.0 && ((*lHisto).second->GetRMS()<15.0)){ 00366 //double t0_estim = (*lHisto).second->GetMean() + 400; 00367 //if(hT0SectorHisto == 0){ 00368 // hT0SectorHisto = new TH1D("hT0AllLayerOfSector","T0 from pulses per layer in sector", 00369 // //20, t0_estim-100, t0_estim+100); 00370 // 700, 0, 7000); 00371 //} 00372 //if(debug) 00373 // cout<<" accepted + 400ns"<<endl; 00374 //hT0SectorHisto->Fill((*lHisto).second->GetMean() + 400); 00375 //} 00376 if(debug) 00377 cout<<endl; 00378 00379 theT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetMean(); 00380 theSigmaT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetRMS();*/ 00381 } 00382 /*if(!hT0SectorHisto){ 00383 cout<<"[DTT0CalibrationNew]: All the t0 per layer are still uncorrect: trying with greater number of events"<<endl; 00384 eventsForLayerT0 = eventsForLayerT0*2; 00385 return; 00386 } 00387 if(debug) 00388 cout<<"[DTT0CalibrationNew] t0 reference for this sector "<< 00389 hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())<<endl;*/ 00390 if(increaseEvtsForLayerT0){ 00391 cout<<"[DTT0CalibrationNew]: t0 per layer are still uncorrect: trying with greater number of events"<<endl; 00392 eventsForLayerT0 = eventsForLayerT0*2; 00393 return; 00394 } 00395 } 00396 }
Compute the mean and the RMS of the t0 from the maps and write them to the DB with channel granularity.
Loop on superlayer to correct between even-odd layers (2 different test pulse lines!)
Change t0 absolute reference -> from sector peak to chamber average
Write the t0 map into DB
Reimplemented from edm::EDAnalyzer.
Definition at line 398 of file DTT0CalibrationNew.cc.
References funct::abs(), DTT0::begin(), DTTimeUnits::counts, GenMuonPlsPt100GeV_cfg::cout, debug, dtGeom, DTT0::end(), lat::endl(), nDigiPerWire, qK, DTT0::set(), sl, funct::sqrt(), theAbsoluteT0PerWire, theCountT0ByChamber, theFile, theHistoLayerMap, theHistoWireMap, theRefT0ByChamber, theRelativeT0PerWire, theSigmaT0PerWire, theSumT0ByChamber, tzero, and DTCalibDBUtils::writeToDB().
00398 { 00399 00400 DTT0* t0s = new DTT0(); 00401 DTT0* t0sWRTChamber = new DTT0(); 00402 00403 if(debug) 00404 cout << "[DTT0CalibrationNewPerLayer]Writing histos to file!" << endl; 00405 00406 theFile->cd(); 00407 //hT0SectorHisto->Write(); 00408 for(map<DTWireId, TH1I*>::const_iterator wHisto = theHistoWireMap.begin(); 00409 wHisto != theHistoWireMap.end(); 00410 wHisto++) { 00411 (*wHisto).second->Write(); 00412 } 00413 for(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin(); 00414 lHisto != theHistoLayerMap.end(); 00415 lHisto++) { 00416 (*lHisto).second->Write(); 00417 } 00418 00419 if(debug) 00420 cout << "[DTT0CalibrationNew] Compute and store t0 and sigma per wire" << endl; 00421 00422 for(map<DTChamberId,double>::const_iterator chamber = theSumT0ByChamber.begin(); 00423 chamber != theSumT0ByChamber.end(); 00424 ++chamber) theRefT0ByChamber[(*chamber).first] = theSumT0ByChamber[(*chamber).first]/((double)theCountT0ByChamber[(*chamber).first]); 00425 00426 for(map<DTWireId, double>::const_iterator wiret0 = theAbsoluteT0PerWire.begin(); 00427 wiret0 != theAbsoluteT0PerWire.end(); 00428 wiret0++){ 00429 if(nDigiPerWire[(*wiret0).first]){ 00430 double t0 = (*wiret0).second/nDigiPerWire[(*wiret0).first]; 00431 DTChamberId chamberId = ((*wiret0).first).chamberId(); 00432 //theRelativeT0PerWire[(*wiret0).first] = t0 - hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin()); 00433 theRelativeT0PerWire[(*wiret0).first] = t0 - theRefT0ByChamber[chamberId]; 00434 cout<<"Wire "<<(*wiret0).first<<" has t0 "<<t0<<"(absolute) "<<theRelativeT0PerWire[(*wiret0).first]<<"(relative)"; 00435 00436 //theSigmaT0PerWire[(*wiret0).first] = sqrt((theSigmaT0PerWire[(*wiret0).first] / nDigiPerWire[(*wiret0).first]) - t0*t0); 00437 theSigmaT0PerWire[(*wiret0).first] = sqrt(qK[(*wiret0).first]/nDigiPerWire[(*wiret0).first]); 00438 cout<<" sigma "<<theSigmaT0PerWire[(*wiret0).first]<<endl; 00439 } 00440 else{ 00441 cout<<"[DTT0CalibrationNew] ERROR: no digis in wire "<<(*wiret0).first<<endl; 00442 abort(); 00443 } 00444 } 00445 00447 // Get all the sls from the setup 00448 const vector<DTSuperLayer*> superLayers = dtGeom->superLayers(); 00449 // Loop over all SLs 00450 for(vector<DTSuperLayer*>::const_iterator sl = superLayers.begin(); 00451 sl != superLayers.end(); sl++) { 00452 00453 00454 //Compute mean for odd and even superlayers 00455 double oddLayersMean=0; 00456 double evenLayersMean=0; 00457 double oddLayersDen=0; 00458 double evenLayersDen=0; 00459 for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin(); 00460 wiret0 != theRelativeT0PerWire.end(); 00461 wiret0++){ 00462 if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){ 00463 if(debug) 00464 cout<<"[DTT0CalibrationNew] Superlayer "<<(*sl)->id() 00465 <<"layer " <<(*wiret0).first.layerId().layer()<<" with "<<(*wiret0).second<<endl; 00466 if(((*wiret0).first.layerId().layer()) % 2){ 00467 oddLayersMean = oddLayersMean + (*wiret0).second; 00468 oddLayersDen++; 00469 } 00470 else{ 00471 evenLayersMean = evenLayersMean + (*wiret0).second; 00472 evenLayersDen++; 00473 } 00474 } 00475 } 00476 oddLayersMean = oddLayersMean/oddLayersDen; 00477 evenLayersMean = evenLayersMean/evenLayersDen; 00478 if(debug && oddLayersMean) 00479 cout<<"[DTT0CalibrationNew] Relative T0 mean for odd layers "<<oddLayersMean<<" even layers"<<evenLayersMean<<endl; 00480 00481 //Compute sigma for odd and even superlayers 00482 double oddLayersSigma=0; 00483 double evenLayersSigma=0; 00484 for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin(); 00485 wiret0 != theRelativeT0PerWire.end(); 00486 wiret0++){ 00487 if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){ 00488 if(((*wiret0).first.layerId().layer()) % 2){ 00489 oddLayersSigma = oddLayersSigma + ((*wiret0).second - oddLayersMean) * ((*wiret0).second - oddLayersMean); 00490 } 00491 else{ 00492 evenLayersSigma = evenLayersSigma + ((*wiret0).second - evenLayersMean) * ((*wiret0).second - evenLayersMean); 00493 } 00494 } 00495 } 00496 oddLayersSigma = oddLayersSigma/oddLayersDen; 00497 evenLayersSigma = evenLayersSigma/evenLayersDen; 00498 oddLayersSigma = sqrt(oddLayersSigma); 00499 evenLayersSigma = sqrt(evenLayersSigma); 00500 00501 if(debug && oddLayersMean) 00502 cout<<"[DTT0CalibrationNew] Relative T0 sigma for odd layers "<<oddLayersSigma<<" even layers"<<evenLayersSigma<<endl; 00503 00504 //Recompute the mean for odd and even superlayers discarding fluctations 00505 double oddLayersFinalMean=0; 00506 double evenLayersFinalMean=0; 00507 for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin(); 00508 wiret0 != theRelativeT0PerWire.end(); 00509 wiret0++){ 00510 if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){ 00511 if(((*wiret0).first.layerId().layer()) % 2){ 00512 if(abs((*wiret0).second - oddLayersMean) < (2*oddLayersSigma)) 00513 oddLayersFinalMean = oddLayersFinalMean + (*wiret0).second; 00514 } 00515 else{ 00516 if(abs((*wiret0).second - evenLayersMean) < (2*evenLayersSigma)) 00517 evenLayersFinalMean = evenLayersFinalMean + (*wiret0).second; 00518 } 00519 } 00520 } 00521 oddLayersFinalMean = oddLayersFinalMean/oddLayersDen; 00522 evenLayersFinalMean = evenLayersFinalMean/evenLayersDen; 00523 if(debug && oddLayersMean) 00524 cout<<"[DTT0CalibrationNew] Final relative T0 mean for odd layers "<<oddLayersFinalMean<<" even layers"<<evenLayersFinalMean<<endl; 00525 00526 for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin(); 00527 wiret0 != theRelativeT0PerWire.end(); 00528 wiret0++){ 00529 if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){ 00530 double t0=-999; 00531 if(((*wiret0).first.layerId().layer()) % 2) 00532 t0 = (*wiret0).second + (evenLayersFinalMean - oddLayersFinalMean); 00533 else 00534 t0 = (*wiret0).second; 00535 00536 cout<<"[DTT0CalibrationNew] Wire "<<(*wiret0).first<<" has t0 "<<(*wiret0).second<<" (relative, after even-odd layer corrections) " 00537 <<" sigma "<<theSigmaT0PerWire[(*wiret0).first]<<endl; 00538 //Store the results into DB 00539 t0s->set((*wiret0).first, t0, theSigmaT0PerWire[(*wiret0).first],DTTimeUnits::counts); 00540 } 00541 } 00542 } 00543 00545 if(debug) 00546 cout << "[DTT0CalibrationNew]Computing relative t0 wrt to chamber average" << endl; 00547 //Compute the reference for each chamber 00548 map<DTChamberId,double> sumT0ByChamber; 00549 map<DTChamberId,int> countT0ByChamber; 00550 for(DTT0::const_iterator tzero = t0s->begin(); 00551 tzero != t0s->end(); tzero++) { 00552 DTChamberId chamberId((*tzero).first.wheelId, 00553 (*tzero).first.stationId, 00554 (*tzero).first.sectorId); 00555 sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + (*tzero).second.t0mean; 00556 countT0ByChamber[chamberId]++; 00557 } 00558 00559 //Change reference for each wire and store the new t0s in the new map 00560 for(DTT0::const_iterator tzero = t0s->begin(); 00561 tzero != t0s->end(); tzero++) { 00562 DTChamberId chamberId((*tzero).first.wheelId, 00563 (*tzero).first.stationId, 00564 (*tzero).first.sectorId); 00565 double t0mean = ((*tzero).second.t0mean) - (sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]); 00566 double t0rms = (*tzero).second.t0rms; 00567 DTWireId wireId((*tzero).first.wheelId, 00568 (*tzero).first.stationId, 00569 (*tzero).first.sectorId, 00570 (*tzero).first.slId, 00571 (*tzero).first.layerId, 00572 (*tzero).first.cellId); 00573 t0sWRTChamber->set(wireId, 00574 t0mean, 00575 t0rms, 00576 DTTimeUnits::counts); 00577 if(debug){ 00578 //cout<<"Chamber "<<chamberId<<" has reference "<<(sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]); 00579 cout<<"Changing t0 of wire "<<wireId<<" from "<<(*tzero).second.t0mean<<" to "<<t0mean<<endl; 00580 } 00581 } 00582 00584 if(debug) 00585 cout << "[DTT0CalibrationNew]Writing values in DB!" << endl; 00586 // FIXME: to be read from cfg? 00587 string t0Record = "DTT0Rcd"; 00588 // Write the t0 map to DB 00589 DTCalibDBUtils::writeToDB(t0Record, t0sWRTChamber); 00590 }
string DTT0CalibrationNew::getHistoName | ( | const DTLayerId & | lId | ) | const [private] |
Definition at line 601 of file DTT0CalibrationNew.cc.
References DTLayerId::layer(), DTChamberId::sector(), DTChamberId::station(), DTSuperLayerId::superlayer(), and DTChamberId::wheel().
00601 { 00602 string histoName; 00603 stringstream theStream; 00604 theStream << "Ch_" << lId.wheel() << "_" << lId.station() << "_" << lId.sector() 00605 << "_SL" << lId.superlayer() << "_L" << lId.layer() <<"_hT0Histo"; 00606 theStream >> histoName; 00607 return histoName; 00608 }
string DTT0CalibrationNew::getHistoName | ( | const DTWireId & | wId | ) | const [private] |
Definition at line 592 of file DTT0CalibrationNew.cc.
References DTLayerId::layer(), DTChamberId::sector(), DTChamberId::station(), DTSuperLayerId::superlayer(), DTChamberId::wheel(), and DTWireId::wire().
Referenced by analyze().
00592 { 00593 string histoName; 00594 stringstream theStream; 00595 theStream << "Ch_" << wId.wheel() << "_" << wId.station() << "_" << wId.sector() 00596 << "_SL" << wId.superlayer() << "_L" << wId.layer() << "_W"<< wId.wire() <<"_hT0Histo"; 00597 theStream >> histoName; 00598 return histoName; 00599 }
std::vector<std::string> DTT0CalibrationNew::cellsWithHistos [private] |
bool DTT0CalibrationNew::debug [private] |
Definition at line 56 of file DTT0CalibrationNew.h.
Referenced by analyze(), DTT0CalibrationNew(), endJob(), and ~DTT0CalibrationNew().
std::string DTT0CalibrationNew::digiLabel [private] |
Definition at line 59 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
edm::ESHandle<DTGeometry> DTT0CalibrationNew::dtGeom [private] |
unsigned int DTT0CalibrationNew::eventsForLayerT0 [private] |
Definition at line 69 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
unsigned int DTT0CalibrationNew::eventsForWireT0 [private] |
Definition at line 71 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
TH1D* DTT0CalibrationNew::hT0SectorHisto [private] |
std::map<DTWireId,double> DTT0CalibrationNew::mK [private] |
std::map<DTWireId,double> DTT0CalibrationNew::mK_ref [private] |
std::map<DTWireId,int> DTT0CalibrationNew::nDigiPerWire [private] |
std::map<DTWireId,int> DTT0CalibrationNew::nDigiPerWire_ref [private] |
unsigned int DTT0CalibrationNew::nevents [private] |
Definition at line 67 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
std::map<DTWireId,double> DTT0CalibrationNew::qK [private] |
unsigned int DTT0CalibrationNew::rejectDigiFromPeak [private] |
Definition at line 83 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
unsigned int DTT0CalibrationNew::retryForLayerT0 [private] |
int DTT0CalibrationNew::selSector [private] |
Definition at line 91 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
int DTT0CalibrationNew::selWheel [private] |
Definition at line 89 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
TSpectrum* DTT0CalibrationNew::spectrum [private] |
Definition at line 98 of file DTT0CalibrationNew.h.
Referenced by analyze(), DTT0CalibrationNew(), and ~DTT0CalibrationNew().
std::map<DTWireId,double> DTT0CalibrationNew::theAbsoluteT0PerWire [private] |
std::string DTT0CalibrationNew::theCalibSector [private] |
Definition at line 90 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
std::string DTT0CalibrationNew::theCalibWheel [private] |
Definition at line 88 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
std::map<DTChamberId,int> DTT0CalibrationNew::theCountT0ByChamber [private] |
TFile* DTT0CalibrationNew::theFile [private] |
Definition at line 62 of file DTT0CalibrationNew.h.
Referenced by analyze(), DTT0CalibrationNew(), endJob(), and ~DTT0CalibrationNew().
std::map<DTLayerId, TH1I*> DTT0CalibrationNew::theHistoLayerMap [private] |
std::map<DTWireId,TH1I*> DTT0CalibrationNew::theHistoWireMap [private] |
TFile* DTT0CalibrationNew::theOutputFile [private] |
Definition at line 64 of file DTT0CalibrationNew.h.
std::map<DTChamberId,double> DTT0CalibrationNew::theRefT0ByChamber [private] |
std::map<DTWireId,double> DTT0CalibrationNew::theRelativeT0PerWire [private] |
std::map<std::string,double> DTT0CalibrationNew::theSigmaT0LayerMap [private] |
Definition at line 117 of file DTT0CalibrationNew.h.
std::map<DTWireId,double> DTT0CalibrationNew::theSigmaT0PerWire [private] |
std::map<DTChamberId,double> DTT0CalibrationNew::theSumT0ByChamber [private] |
std::map<std::string,double> DTT0CalibrationNew::theT0LayerMap [private] |
Definition at line 116 of file DTT0CalibrationNew.h.
std::map<DTLayerId,double> DTT0CalibrationNew::theTPPeakMap [private] |
unsigned int DTT0CalibrationNew::timeBoxWidth [private] |
Definition at line 80 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
double DTT0CalibrationNew::tpPeakWidth [private] |
Definition at line 74 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
double DTT0CalibrationNew::tpPeakWidthPerLayer [private] |
Definition at line 77 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().
std::vector<DTWireId> DTT0CalibrationNew::wireIdWithHistos [private] |
Definition at line 101 of file DTT0CalibrationNew.h.
Referenced by analyze(), and DTT0CalibrationNew().