26 #include "TSpectrum.h"
38 cout <<
"[DTT0CalibrationNew]Constructor called!" << endl;
47 theFile =
new TFile(rootFileName.c_str(),
"RECREATE");
50 if(theCalibWheel !=
"All") {
53 linestr << theCalibWheel;
55 cout <<
"[DTT0CalibrationNewPerLayer] chosen wheel " << selWheel << endl;
60 if(theCalibSector !=
"All") {
63 linestr << theCalibSector;
65 cout <<
"[DTT0CalibrationNewPerLayer] chosen sector " << selSector << endl;
68 vector<string> defaultCell;
69 defaultCell.push_back(
"None");
71 for(vector<string>::const_iterator cell = cellsWithHistos.begin(); cell != cellsWithHistos.end(); cell++){
72 if((*cell) !=
"None"){
74 int wheel,sector,
station,sl,layer,wire;
76 linestr >> wheel >> sector >> station >> sl >> layer >> wire;
77 wireIdWithHistos.push_back(
DTWireId(wheel,station,sector,sl,layer,wire));
84 eventsForLayerT0 = pset.
getParameter<
unsigned int>(
"eventsForLayerT0");
85 eventsForWireT0 = pset.
getParameter<
unsigned int>(
"eventsForWireT0");
87 tpPeakWidthPerLayer = pset.
getParameter<
double>(
"tpPeakWidthPerLayer");
88 timeBoxWidth = pset.
getParameter<
unsigned int>(
"timeBoxWidth");
89 rejectDigiFromPeak = pset.
getParameter<
unsigned int>(
"rejectDigiFromPeak");
91 spectrum =
new TSpectrum(5);
98 cout <<
"[DTT0CalibrationNew]Destructor called!" << endl;
107 cout <<
"--- [DTT0CalibrationNew] Analysing Event: #Run: " << event.
id().
run()
108 <<
" #Event: " <<
event.id().event() << endl;
124 for (dtLayerIt = digis->begin();
125 dtLayerIt != digis->end();
131 const DTLayerId layerId = (*dtLayerIt).first;
143 float tTrig,tTrigRMS, kFactor;
147 <<
" tTrig,tTrigRMS= " << tTrig <<
", " << tTrigRMS << endl;
152 digi != digiRange.second;
154 double t0 = (*digi).countsTDC();
157 if(
nevents < eventsForLayerT0){
159 TH1I *hT0LayerHisto = theHistoLayerMap[layerId];
161 if(hT0LayerHisto == 0){
163 float hT0Min = tTrig - 2*tTrigRMS;
164 float hT0Max = hT0Min + timeBoxWidth;
166 "T0 from pulses by layer (TDC counts, 1 TDC count = 0.781 ns)",
167 timeBoxWidth,hT0Min,hT0Max);
169 cout <<
" New T0 per Layer Histo: " << hT0LayerHisto->GetName() << endl;
170 theHistoLayerMap[layerId] = hT0LayerHisto;
175 if(hT0LayerHisto != 0) {
178 hT0LayerHisto->Fill(t0);
185 const DTWireId wireId(layerId, (*digi).wire());
187 cout <<
" Wire: " << wireId << endl
188 <<
" time (TDC counts): " << (*digi).countsTDC()<< endl;
192 vector<DTWireId>::iterator it =
find(wireIdWithHistos.begin(),wireIdWithHistos.end(),wireId);
193 if (it!=wireIdWithHistos.end()){
195 TH1I *hT0WireHisto = theHistoWireMap[wireId];
197 if(hT0WireHisto == 0){
199 hT0WireHisto =
new TH1I(
getHistoName(wireId).c_str(),
"T0 from pulses by wire (TDC counts, 1 TDC count = 0.781 ns)",7000,0,7000);
203 cout <<
" New T0 per wire Histo: " << hT0WireHisto->GetName() << endl;
204 theHistoWireMap[wireId] = hT0WireHisto;
208 if(hT0WireHisto != 0) {
211 hT0WireHisto->Fill(t0);
229 if(fabs(theTPPeakMap[layerId] - t0) > rejectDigiFromPeak){
231 cout<<
"digi skipped because t0 too far from peak " << theTPPeakMap[layerId] << endl;
236 theSumT0ByChamber[chamberId] = theSumT0ByChamber[chamberId] + t0;
237 theCountT0ByChamber[chamberId]++;
240 if(
nevents< (eventsForLayerT0 + eventsForWireT0)){
241 if(!nDigiPerWire_ref[wireId]){
244 nDigiPerWire_ref[wireId] = nDigiPerWire_ref[wireId] + 1;
245 mK_ref[wireId] = mK_ref[wireId] + (t0-mK_ref[wireId])/nDigiPerWire_ref[wireId];
248 else if(
nevents>(eventsForLayerT0 + eventsForWireT0)){
249 if(
abs(t0-mK_ref[wireId]) > tpPeakWidth)
251 if(!nDigiPerWire[wireId]){
252 theAbsoluteT0PerWire[wireId] = 0;
256 nDigiPerWire[wireId] = nDigiPerWire[wireId] + 1;
257 theAbsoluteT0PerWire[wireId] = theAbsoluteT0PerWire[wireId] + t0;
259 qK[wireId] = qK[wireId] + ((nDigiPerWire[wireId]-1)*(t0-mK[wireId])*(t0-mK[wireId])/nDigiPerWire[wireId]);
260 mK[wireId] = mK[wireId] + (t0-mK[wireId])/nDigiPerWire[wireId];
267 if(
nevents == eventsForLayerT0){
268 bool increaseEvtsForLayerT0 =
false;
269 for(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin();
270 lHisto != theHistoLayerMap.end();
273 cout<<
"Reading histogram "<<(*lHisto).second->GetName()<<
" with mean "<<(*lHisto).second->GetMean()<<
" and RMS "<<(*lHisto).second->GetRMS() << endl;
278 int npeaks = spectrum->Search((*lHisto).second,(tpPeakWidthPerLayer/2.),
"",0.3);
280 float *peaks = spectrum->GetPositionX();
282 vector<float> peakMeans(peaks,peaks + npeaks);
284 sort(peakMeans.begin(),peakMeans.end());
287 float tTrig,tTrigRMS, kFactor;
288 tTrigMap->get((*lHisto).first.superlayerId(), tTrig, tTrigRMS, kFactor,
DTTimeUnits::counts );
290 float timeBoxCenter = (2*tTrig + (float)timeBoxWidth)/2.;
291 float hMin = (*lHisto).second->GetXaxis()->GetXmin();
292 float hMax = (*lHisto).second->GetXaxis()->GetXmax();
293 vector<float>::const_iterator tpPeak = peakMeans.end();
294 for(vector<float>::const_iterator it = peakMeans.begin(); it != peakMeans.end(); ++it){
297 int bin = (*lHisto).second->GetXaxis()->FindBin(mean);
298 float yp = (*lHisto).second->GetBinContent(bin);
299 if(
debug)
cout <<
"Peak : (" << mean <<
"," << yp <<
")" << endl;
302 float previous_peak = (it == peakMeans.begin())?hMin:*(it - 1);
303 float next_peak = (it == (peakMeans.end()-1))?hMax:*(it + 1);
305 float rangemin = mean - (mean - previous_peak)/8.;
306 float rangemax = mean + (next_peak -
mean)/8.;
307 int binmin = (*lHisto).second->GetXaxis()->FindBin(rangemin);
308 int binmax = (*lHisto).second->GetXaxis()->FindBin(rangemax);
309 (*lHisto).second->GetXaxis()->SetRange(binmin,binmax);
311 float rms_seed = (*lHisto).second->GetRMS();
332 if(rms_seed > tpPeakWidthPerLayer)
continue;
334 if(fabs(mean - timeBoxCenter) < fabs(*tpPeak - timeBoxCenter)) tpPeak = it;
345 float selPeak = (tpPeak != peakMeans.end())?*tpPeak:(*lHisto).second->GetBinCenter((*lHisto).second->GetMaximumBin());
346 if(
debug)
cout <<
"Peak selected at " << selPeak << endl;
348 theTPPeakMap[(*lHisto).first] = selPeak;
391 if(increaseEvtsForLayerT0){
392 cout<<
"[DTT0CalibrationNew]: t0 per layer are still uncorrect: trying with greater number of events"<<endl;
393 eventsForLayerT0 = eventsForLayerT0*2;
405 cout <<
"[DTT0CalibrationNewPerLayer]Writing histos to file!" << endl;
409 for(map<DTWireId, TH1I*>::const_iterator wHisto = theHistoWireMap.begin();
410 wHisto != theHistoWireMap.end();
412 (*wHisto).second->Write();
414 for(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin();
415 lHisto != theHistoLayerMap.end();
417 (*lHisto).second->Write();
421 cout <<
"[DTT0CalibrationNew] Compute and store t0 and sigma per wire" << endl;
423 for(map<DTChamberId,double>::const_iterator chamber = theSumT0ByChamber.begin();
424 chamber != theSumT0ByChamber.end();
425 ++chamber) theRefT0ByChamber[(*chamber).first] = theSumT0ByChamber[(*chamber).first]/((
double)theCountT0ByChamber[(*chamber).first]);
427 for(map<DTWireId, double>::const_iterator wiret0 = theAbsoluteT0PerWire.begin();
428 wiret0 != theAbsoluteT0PerWire.end();
430 if(nDigiPerWire[(*wiret0).first]){
431 double t0 = (*wiret0).second/nDigiPerWire[(*wiret0).first];
432 DTChamberId chamberId = ((*wiret0).first).chamberId();
434 theRelativeT0PerWire[(*wiret0).first] = t0 - theRefT0ByChamber[chamberId];
435 cout<<
"Wire "<<(*wiret0).first<<
" has t0 "<<t0<<
"(absolute) "<<theRelativeT0PerWire[(*wiret0).first]<<
"(relative)";
438 theSigmaT0PerWire[(*wiret0).first] =
sqrt(qK[(*wiret0).first]/nDigiPerWire[(*wiret0).first]);
439 cout<<
" sigma "<<theSigmaT0PerWire[(*wiret0).first]<<endl;
442 cout<<
"[DTT0CalibrationNew] ERROR: no digis in wire "<<(*wiret0).first<<endl;
449 const vector<DTSuperLayer*> superLayers = dtGeom->superLayers();
451 for(vector<DTSuperLayer*>::const_iterator sl = superLayers.begin();
452 sl != superLayers.end(); sl++) {
456 double oddLayersMean=0;
457 double evenLayersMean=0;
458 double oddLayersDen=0;
459 double evenLayersDen=0;
460 for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
461 wiret0 != theRelativeT0PerWire.end();
463 if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
465 cout<<
"[DTT0CalibrationNew] Superlayer "<<(*sl)->id()
466 <<
"layer " <<(*wiret0).first.layerId().layer()<<
" with "<<(*wiret0).second<<endl;
467 if(((*wiret0).first.layerId().layer()) % 2){
468 oddLayersMean = oddLayersMean + (*wiret0).second;
472 evenLayersMean = evenLayersMean + (*wiret0).second;
477 oddLayersMean = oddLayersMean/oddLayersDen;
478 evenLayersMean = evenLayersMean/evenLayersDen;
479 if(
debug && oddLayersMean)
480 cout<<
"[DTT0CalibrationNew] Relative T0 mean for odd layers "<<oddLayersMean<<
" even layers"<<evenLayersMean<<endl;
483 double oddLayersSigma=0;
484 double evenLayersSigma=0;
485 for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
486 wiret0 != theRelativeT0PerWire.end();
488 if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
489 if(((*wiret0).first.layerId().layer()) % 2){
490 oddLayersSigma = oddLayersSigma + ((*wiret0).second - oddLayersMean) * ((*wiret0).second - oddLayersMean);
493 evenLayersSigma = evenLayersSigma + ((*wiret0).second - evenLayersMean) * ((*wiret0).second - evenLayersMean);
497 oddLayersSigma = oddLayersSigma/oddLayersDen;
498 evenLayersSigma = evenLayersSigma/evenLayersDen;
499 oddLayersSigma =
sqrt(oddLayersSigma);
500 evenLayersSigma =
sqrt(evenLayersSigma);
502 if(
debug && oddLayersMean)
503 cout<<
"[DTT0CalibrationNew] Relative T0 sigma for odd layers "<<oddLayersSigma<<
" even layers"<<evenLayersSigma<<endl;
506 double oddLayersFinalMean=0;
507 double evenLayersFinalMean=0;
508 for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
509 wiret0 != theRelativeT0PerWire.end();
511 if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
512 if(((*wiret0).first.layerId().layer()) % 2){
513 if(
abs((*wiret0).second - oddLayersMean) < (2*oddLayersSigma))
514 oddLayersFinalMean = oddLayersFinalMean + (*wiret0).second;
517 if(
abs((*wiret0).second - evenLayersMean) < (2*evenLayersSigma))
518 evenLayersFinalMean = evenLayersFinalMean + (*wiret0).second;
522 oddLayersFinalMean = oddLayersFinalMean/oddLayersDen;
523 evenLayersFinalMean = evenLayersFinalMean/evenLayersDen;
524 if(
debug && oddLayersMean)
525 cout<<
"[DTT0CalibrationNew] Final relative T0 mean for odd layers "<<oddLayersFinalMean<<
" even layers"<<evenLayersFinalMean<<endl;
527 for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
528 wiret0 != theRelativeT0PerWire.end();
530 if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
532 if(((*wiret0).first.layerId().layer()) % 2)
533 t0 = (*wiret0).second + (evenLayersFinalMean - oddLayersFinalMean);
535 t0 = (*wiret0).second;
537 cout<<
"[DTT0CalibrationNew] Wire "<<(*wiret0).first<<
" has t0 "<<(*wiret0).second<<
" (relative, after even-odd layer corrections) "
538 <<
" sigma "<<theSigmaT0PerWire[(*wiret0).first]<<endl;
547 cout <<
"[DTT0CalibrationNew]Computing relative t0 wrt to chamber average" << endl;
549 map<DTChamberId,double> sumT0ByChamber;
550 map<DTChamberId,int> countT0ByChamber;
558 int channelId =
tzero->channelId;
559 if ( channelId == 0 )
continue;
567 sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + t0mean_f;
569 countT0ByChamber[chamberId]++;
587 int channelId =
tzero->channelId;
588 if ( channelId == 0 )
continue;
597 double t0mean = t0mean_f - (sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
598 double t0rms = t0rms_f;
600 t0sWRTChamber->
set(wireId,
607 cout<<
"Changing t0 of wire "<<wireId<<
" from "<<
tzero->t0mean<<
" to "<<t0mean<<endl;
613 cout <<
"[DTT0CalibrationNew]Writing values in DB!" << endl;
615 string t0Record =
"DTT0Rcd";
622 stringstream theStream;
625 theStream >> histoName;
631 stringstream theStream;
634 theStream >> histoName;
const_iterator begin() const
int set(int wheelId, int stationId, int sectorId, int slId, int layerId, int cellId, float t0mean, float t0rms, DTTimeUnits::type unit)
T getParameter(std::string const &) const
EventNumber_t event() const
T getUntrackedParameter(std::string const &, T const &) const
DTChamberId chamberId() const
Return the corresponding ChamberId.
void endJob()
Compute the mean and the RMS of the t0 from the maps and write them to the DB with channel granularit...
int layer() const
Return the layer number.
DTSuperLayerId superlayerId() const
Return the corresponding SuperLayerId.
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
int get(int wheelId, int stationId, int sectorId, int slId, int layerId, int cellId, float &t0mean, float &t0rms, DTTimeUnits::type unit) const
std::vector< DTT0Data >::const_iterator const_iterator
Access methods to data.
virtual ~DTT0CalibrationNew()
Destructor.
Abs< T >::type abs(const T &t)
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
int wire() const
Return the wire number.
const_iterator end() const
int superlayer() const
Return the superlayer number (deprecated method name)
DTT0CalibrationNew(const edm::ParameterSet &pset)
Constructor.
std::vector< DTDigi >::const_iterator const_iterator
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup)
Fill the maps with t0 (by channel)
static const double tzero[3]
std::pair< const_iterator, const_iterator > Range
std::string getHistoName(const DTWireId &wId) const
int station() const
Return the station number.
int wheel() const
Return the wheel number.
static void writeToDB(std::string record, T *payload)