CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Private Attributes
DTT0CalibrationNew Class Reference

#include <DTT0CalibrationNew.h>

Inheritance diagram for DTT0CalibrationNew:
edm::EDAnalyzer edm::EDConsumerBase

Public Member Functions

void analyze (const edm::Event &event, const edm::EventSetup &eventSetup)
 Fill the maps with t0 (by channel) More...
 
 DTT0CalibrationNew (const edm::ParameterSet &pset)
 Constructor. More...
 
void endJob ()
 Compute the mean and the RMS of the t0 from the maps and write them to the DB with channel granularity. More...
 
virtual ~DTT0CalibrationNew ()
 Destructor. More...
 
- Public Member Functions inherited from edm::EDAnalyzer
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
 EDAnalyzer ()
 
ModuleDescription const & moduleDescription () const
 
std::string workerType () const
 
virtual ~EDAnalyzer ()
 
- Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfoconsumesInfo () const
 
 EDConsumerBase ()
 
 EDConsumerBase (EDConsumerBase const &)=delete
 
 EDConsumerBase (EDConsumerBase &&)=default
 
ProductResolverIndexAndSkipBit indexFrom (EDGetToken, BranchType, TypeID const &) const
 
void itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
void itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
std::vector< ProductResolverIndexAndSkipBit > const & itemsToGetFrom (BranchType iType) const
 
void labelsForToken (EDGetToken iToken, Labels &oLabels) const
 
void modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const
 
EDConsumerBase const & operator= (EDConsumerBase const &)=delete
 
EDConsumerBaseoperator= (EDConsumerBase &&)=default
 
bool registeredToConsume (ProductResolverIndex, bool, BranchType) const
 
bool registeredToConsumeMany (TypeID const &, BranchType) const
 
void updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)
 
virtual ~EDConsumerBase () noexcept(false)
 

Private Member Functions

std::string getHistoName (const DTWireId &wId) const
 
std::string getHistoName (const DTLayerId &lId) const
 

Private Attributes

std::vector< std::string > cellsWithHistos
 
std::string dbLabel
 
bool debug
 
std::string digiLabel
 
edm::ESHandle< DTGeometrydtGeom
 
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< DTWireIdwireIdWithHistos
 

Additional Inherited Members

- Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer ModuleType
 
- Public Types inherited from edm::EDConsumerBase
typedef ProductLabels Labels
 
- Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &)
 
- Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > consumes (edm::InputTag const &tag)
 
EDGetToken consumes (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken consumes (TypeToGet const &id, edm::InputTag const &tag)
 
ConsumesCollector consumesCollector ()
 Use a ConsumesCollector to gather consumes information from helper functions. More...
 
template<typename ProductType , BranchType B = InEvent>
void consumesMany ()
 
void consumesMany (const TypeToGet &id)
 
template<BranchType B>
void consumesMany (const TypeToGet &id)
 
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > mayConsume (edm::InputTag const &tag)
 
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 

Detailed Description

Analyzer class computes the mean and RMS of t0 from pulses. 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

Author
S. Bolognesi - INFN Torino

Definition at line 29 of file DTT0CalibrationNew.h.

Constructor & Destructor Documentation

DTT0CalibrationNew::DTT0CalibrationNew ( const edm::ParameterSet pset)

Constructor.

Definition at line 34 of file DTT0CalibrationNew.cc.

References gather_cfg::cout, RecoTauHPSTancTauProdcuer_cfi::dbLabel, debug, CastorSimpleReconstructor_cfi::digiLabel, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), DTAnalyzerDetailed_cfi::rootFileName, relativeConstraints::station, interactiveExample::theFile, and makeMuonMisalignmentScenario::wheel.

34  {
35  // Get the debug parameter for verbose output
36  debug = pset.getUntrackedParameter<bool>("debug");
37  if(debug)
38  cout << "[DTT0CalibrationNew]Constructor called!" << endl;
39 
40  // Get the label to retrieve digis from the event
41  digiLabel = pset.getUntrackedParameter<string>("digiLabel");
42 
43  dbLabel = pset.getUntrackedParameter<string>("dbLabel", "");
44 
45  // The root file which contain the histos per layer
46  string rootFileName = pset.getUntrackedParameter<string>("rootFileName","DTT0PerLayer.root");
47  theFile = new TFile(rootFileName.c_str(), "RECREATE");
48 
49  theCalibWheel = pset.getUntrackedParameter<string>("calibWheel", "All"); //FIXME amke a vector of integer instead of a string
50  if(theCalibWheel != "All") {
51  stringstream linestr;
52  int selWheel;
53  linestr << theCalibWheel;
54  linestr >> selWheel;
55  cout << "[DTT0CalibrationNewPerLayer] chosen wheel " << selWheel << endl;
56  }
57 
58  // Sector/s to calibrate
59  theCalibSector = pset.getUntrackedParameter<string>("calibSector", "All"); //FIXME amke a vector of integer instead of a string
60  if(theCalibSector != "All") {
61  stringstream linestr;
62  int selSector;
63  linestr << theCalibSector;
64  linestr >> selSector;
65  cout << "[DTT0CalibrationNewPerLayer] chosen sector " << selSector << endl;
66  }
67 
68  vector<string> defaultCell;
69  defaultCell.push_back("None");
70  cellsWithHistos = pset.getUntrackedParameter<vector<string> >("cellsWithHisto", defaultCell);
71  for(vector<string>::const_iterator cell = cellsWithHistos.begin(); cell != cellsWithHistos.end(); ++cell){
72  if((*cell) != "None"){
73  stringstream linestr;
74  int wheel,sector,station,sl,layer,wire;
75  linestr << (*cell);
76  linestr >> wheel >> sector >> station >> sl >> layer >> wire;
77  wireIdWithHistos.push_back(DTWireId(wheel,station,sector,sl,layer,wire));
78  }
79  }
80 
82 
83  nevents=0;
84  eventsForLayerT0 = pset.getParameter<unsigned int>("eventsForLayerT0");
85  eventsForWireT0 = pset.getParameter<unsigned int>("eventsForWireT0");
86  tpPeakWidth = pset.getParameter<double>("tpPeakWidth");
87  tpPeakWidthPerLayer = pset.getParameter<double>("tpPeakWidthPerLayer");
88  timeBoxWidth = pset.getParameter<unsigned int>("timeBoxWidth");
89  rejectDigiFromPeak = pset.getParameter<unsigned int>("rejectDigiFromPeak");
90 
91  spectrum = new TSpectrum(5);
92  retryForLayerT0 = 0;
93 }
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
unsigned int eventsForLayerT0
unsigned int retryForLayerT0
unsigned int eventsForWireT0
unsigned int timeBoxWidth
std::vector< std::string > cellsWithHistos
std::vector< DTWireId > wireIdWithHistos
unsigned int rejectDigiFromPeak
std::string theCalibSector
DTT0CalibrationNew::~DTT0CalibrationNew ( )
virtual

Destructor.

Definition at line 96 of file DTT0CalibrationNew.cc.

References gather_cfg::cout, debug, and interactiveExample::theFile.

96  {
97  if(debug)
98  cout << "[DTT0CalibrationNew]Destructor called!" << endl;
99 
100  delete spectrum;
101  theFile->Close();
102 }

Member Function Documentation

void DTT0CalibrationNew::analyze ( const edm::Event event,
const edm::EventSetup eventSetup 
)

Fill the maps with t0 (by channel)

Perform the real analysis.

Definition at line 105 of file DTT0CalibrationNew.cc.

References funct::abs(), stringResolutionProvider_cfi::bin, DTSuperLayerId::chamberId(), DTTimeUnits::counts, gather_cfg::cout, RecoTauHPSTancTauProdcuer_cfi::dbLabel, debug, CastorSimpleReconstructor_cfi::digiLabel, edm::EventID::event(), spr::find(), objects.autophobj::float, edm::EventSetup::get(), DTTtrig::get(), mergeVDriftHistosByStation::getHistoName(), edm::EventBase::id(), DTAnalyzerDetailed_cfi::kFactor, RecoTauDiscriminantConfiguration::mean, edm::EventID::run(), DTChamberId::sector(), DTLayerId::superlayerId(), cscNeutronWriter_cfi::t0, interactiveExample::theFile, and DTChamberId::wheel().

105  {
106  if(debug || event.id().event() % 500==0)
107  cout << "--- [DTT0CalibrationNew] Analysing Event: #Run: " << event.id().run()
108  << " #Event: " << event.id().event() << endl;
109  nevents++;
110 
111  // Get the digis from the event
113  event.getByLabel(digiLabel, digis);
114 
115  // Get the DT Geometry
116  eventSetup.get<MuonGeometryRecord>().get(dtGeom);
117 
118  // Get ttrig DB
119  edm::ESHandle<DTTtrig> tTrigMap;
120  eventSetup.get<DTTtrigRcd>().get(dbLabel,tTrigMap);
121 
122  // Iterate through all digi collections ordered by LayerId
124  for (dtLayerIt = digis->begin();
125  dtLayerIt != digis->end();
126  ++dtLayerIt){
127  // Get the iterators over the digis associated with this LayerId
128  const DTDigiCollection::Range& digiRange = (*dtLayerIt).second;
129 
130  // Get the layerId
131  const DTLayerId layerId = (*dtLayerIt).first; //FIXME: check to be in the right sector
132  const DTChamberId chamberId = layerId.superlayerId().chamberId();
133 
134  if((theCalibWheel != "All") && (layerId.superlayerId().chamberId().wheel() != selWheel))
135  continue;
136  if((theCalibSector != "All") && (layerId.superlayerId().chamberId().sector() != selSector))
137  continue;
138 
139  //if(debug) {
140  // cout << "Layer " << layerId<<" with "<<distance(digiRange.first, digiRange.second)<<" digi"<<endl;
141  //}
142 
143  float tTrig,tTrigRMS, kFactor;
144  tTrigMap->get(layerId.superlayerId(), tTrig, tTrigRMS, kFactor, DTTimeUnits::counts );
145  if(debug&&(nevents <= 1)){
146  cout << " Superlayer: " << layerId.superlayerId() << endl
147  << " tTrig,tTrigRMS= " << tTrig << ", " << tTrigRMS << endl;
148  }
149 
150  // Loop over all digis in the given layer
151  for (DTDigiCollection::const_iterator digi = digiRange.first;
152  digi != digiRange.second;
153  ++digi) {
154  double t0 = (*digi).countsTDC();
155 
156  //Use first bunch of events to fill t0 per layer
158  //Get the per-layer histo from the map
159  TH1I *hT0LayerHisto = theHistoLayerMap[layerId];
160  //If it doesn't exist, book it
161  if(hT0LayerHisto == 0){
162  theFile->cd();
163  float hT0Min = tTrig - 2*tTrigRMS;
164  float hT0Max = hT0Min + timeBoxWidth;
165  hT0LayerHisto = new TH1I(getHistoName(layerId).c_str(),
166  "T0 from pulses by layer (TDC counts, 1 TDC count = 0.781 ns)",
167  timeBoxWidth,hT0Min,hT0Max);
168  if(debug)
169  cout << " New T0 per Layer Histo: " << hT0LayerHisto->GetName() << endl;
170  theHistoLayerMap[layerId] = hT0LayerHisto;
171  }
172 
173  //Fill the histos
174  theFile->cd();
175  if(hT0LayerHisto != 0) {
176  // if(debug)
177  // cout<<"Filling histo "<<hT0LayerHisto->GetName()<<" with digi "<<t0<<" TDC counts"<<endl;
178  hT0LayerHisto->Fill(t0);
179  }
180  }
181 
182  //Use all the remaining events to compute t0 per wire
184  // Get the wireId
185  const DTWireId wireId(layerId, (*digi).wire());
186  if(debug) {
187  cout << " Wire: " << wireId << endl
188  << " time (TDC counts): " << (*digi).countsTDC()<< endl;
189  }
190 
191  //Fill the histos per wire for the chosen cells
192  vector<DTWireId>::iterator it = find(wireIdWithHistos.begin(),wireIdWithHistos.end(),wireId);
193  if (it!=wireIdWithHistos.end()){
194  //Get the per-wire histo from the map
195  TH1I *hT0WireHisto = theHistoWireMap[wireId];
196  //If it doesn't exist, book it
197  if(hT0WireHisto == 0){
198  theFile->cd();
199  hT0WireHisto = new TH1I(getHistoName(wireId).c_str(),"T0 from pulses by wire (TDC counts, 1 TDC count = 0.781 ns)",7000,0,7000);
200  //hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())-100,
201  //hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())+100);
202  if(debug)
203  cout << " New T0 per wire Histo: " << hT0WireHisto->GetName() << endl;
204  theHistoWireMap[wireId] = hT0WireHisto;
205  }
206  //Fill the histos
207  theFile->cd();
208  if(hT0WireHisto != 0) {
209  //if(debug)
210  // cout<<"Filling histo "<<hT0WireHisto->GetName()<<" with digi "<<t0<<" TDC counts"<<endl;
211  hT0WireHisto->Fill(t0);
212  }
213  }
214 
215  //Check the tzero has reasonable value
216  //float hT0Min = tTrig - 2*tTrigRMS;
217  //float hT0Max = hT0Min + timeBoxWidth;
218  /*if(abs(hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin()) - t0) > rejectDigiFromPeak){
219  if(debug)
220  cout<<"digi skipped because t0 per sector "<<hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())<<endl;
221  continue;
222  }*/
223  /*if((t0 < hT0Min)||(t0 > hT0Max)){
224  if(debug)
225  cout<<"digi skipped because t0 outside of interval (" << hT0Min << "," << hT0Max << ")" <<endl;
226  continue;
227  }*/
228  //Select per layer
229  if(fabs(theTPPeakMap[layerId] - t0) > rejectDigiFromPeak){
230  if(debug)
231  cout<<"digi skipped because t0 too far from peak " << theTPPeakMap[layerId] << endl;
232  continue;
233  }
234 
235  //Find to ref. per chamber
236  theSumT0ByChamber[chamberId] = theSumT0ByChamber[chamberId] + t0;
237  theCountT0ByChamber[chamberId]++;
238 
239  //Use second bunch of events to compute a t0 reference per wire
241  if(!nDigiPerWire_ref[wireId]){
242  mK_ref[wireId] = 0;
243  }
244  nDigiPerWire_ref[wireId] = nDigiPerWire_ref[wireId] + 1;
245  mK_ref[wireId] = mK_ref[wireId] + (t0-mK_ref[wireId])/nDigiPerWire_ref[wireId];
246  }
247  //Use last all the remaining events to compute the mean and sigma t0 per wire
249  if(abs(t0-mK_ref[wireId]) > tpPeakWidth)
250  continue;
251  if(!nDigiPerWire[wireId]){
252  theAbsoluteT0PerWire[wireId] = 0;
253  qK[wireId] = 0;
254  mK[wireId] = 0;
255  }
256  nDigiPerWire[wireId] = nDigiPerWire[wireId] + 1;
257  theAbsoluteT0PerWire[wireId] = theAbsoluteT0PerWire[wireId] + t0;
258  //theSigmaT0PerWire[wireId] = theSigmaT0PerWire[wireId] + (t0*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];
261  }
262  }//end if(nevents>1000)
263  }//end loop on digi
264  }//end loop on layer
265 
266  //Use the t0 per layer histos to have an indication about the t0 position
267  if(nevents == eventsForLayerT0){
268  bool increaseEvtsForLayerT0 = false;
269  for(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin();
270  lHisto != theHistoLayerMap.end();
271  ++lHisto){
272  if(debug)
273  cout<<"Reading histogram "<<(*lHisto).second->GetName()<<" with mean "<<(*lHisto).second->GetMean()<<" and RMS "<<(*lHisto).second->GetRMS() << endl;
274 
275  //Find peaks
276  //int npeaks = spectrum->Search((*lHisto).second,0.5,"goff");
277  //int npeaks = spectrum->Search((*lHisto).second,(tpPeakWidthPerLayer/2.),"goff",0.3);
278  int npeaks = spectrum->Search((*lHisto).second,(tpPeakWidthPerLayer/2.),"",0.3);
279 
280  double *peaks = spectrum->GetPositionX();
281  //Put in a std::vector<float>
282  vector<float> peakMeans(peaks,peaks + npeaks);
283  //Sort the peaks in ascending order
284  sort(peakMeans.begin(),peakMeans.end());
285 
286  //Find best peak -- preliminary criteria: find peak closest to center of time box
287  float tTrig,tTrigRMS, kFactor;
288  tTrigMap->get((*lHisto).first.superlayerId(), tTrig, tTrigRMS, kFactor, DTTimeUnits::counts );
289 
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){
295  float mean = *it;
296 
297  int bin = (*lHisto).second->GetXaxis()->FindBin(mean);
298  float yp = (*lHisto).second->GetBinContent(bin);
299  if(debug) cout << "Peak : (" << mean << "," << yp << ")" << endl;
300 
301  //Find RMS
302  float previous_peak = (it == peakMeans.begin())?hMin:*(it - 1);
303  float next_peak = (it == (peakMeans.end()-1))?hMax:*(it + 1);
304 
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);
310  //RMS estimate
311  float rms_seed = (*lHisto).second->GetRMS();
312 
313  /*rangemin = mean - 2*rms_seed;
314  rangemax = mean + 2*rms_seed;
315  if(debug) cout << "Seed for RMS, Fit min, Fit max: " << rms_seed << ", " << rangemin << ", " << rangemax << endl;
316  //Fit to gaussian
317  string funcname("fitFcn_");
318  funcname += (*lHisto).second->GetName();
319  if(debug) cout << "Fitting function " << funcname << endl;
320  TF1* func = new TF1(funcname.c_str(),"gaus",rangemin,rangemax);
321  func->SetParameters(yp,mean,rms_seed);
322  (*lHisto).second->Fit(func,"Q","",rangemin,rangemax);
323 
324  float fitconst = func->GetParameter(0);
325  float fitmean = func->GetParameter(1);
326  float fitrms = func->GetParameter(2);
327  float chisquare = func->GetChisquare()/func->GetNDF();
328  if(debug) cout << "Gaussian fit constant,mean,RMS,chi2= " << fitconst << ", " << fitmean << ", " << fitrms << ", " << chisquare << endl;*/
329 
330  //Reject peaks with RMS larger than specified
331  //if(fitrms > tpPeakWidth) continue;
332  if(rms_seed > tpPeakWidthPerLayer) continue;
333 
334  if(fabs(mean - timeBoxCenter) < fabs(*tpPeak - timeBoxCenter)) tpPeak = it;
335  }
336  //Didn't find peak
337  /*if(tpPeak == peakMeans.end()){
338  if(retryForLayerT0 < 2){
339  increaseEvtsForLayerT0 = true;
340  retryForLayerT0++;
341  break;
342  }
343  }*/
344 
345  float selPeak = (tpPeak != peakMeans.end())?*tpPeak:(*lHisto).second->GetBinCenter((*lHisto).second->GetMaximumBin());
346  if(debug) cout << "Peak selected at " << selPeak << endl;
347 
348  theTPPeakMap[(*lHisto).first] = selPeak;
349 
350  //Take the mean as a first t0 estimation
351  /*if((*lHisto).second->GetRMS() < tpPeakWidth){
352  if(hT0SectorHisto == 0){
353  hT0SectorHisto = new TH1D("hT0AllLayerOfSector","T0 from pulses per layer in sector",
354  //20, (*lHisto).second->GetMean()-100, (*lHisto).second->GetMean()+100);
355  700, 0, 7000);
356  //300,3300,3600);
357  }
358  if(debug)
359  cout<<" accepted"<<endl;
360  //TH1I* aux_histo = (*lHisto).second;
361  //aux_histo->GetXaxis()->SetRangeUser(3300,3600);
362  hT0SectorHisto->Fill((*lHisto).second->GetMean());
363  //hT0SectorHisto->Fill(aux_histo->GetMean());
364  }
365  //Take the mean of noise + 400ns as a first t0 estimation
366  //if((*lHisto).second->GetRMS()>10.0 && ((*lHisto).second->GetRMS()<15.0)){
367  //double t0_estim = (*lHisto).second->GetMean() + 400;
368  //if(hT0SectorHisto == 0){
369  // hT0SectorHisto = new TH1D("hT0AllLayerOfSector","T0 from pulses per layer in sector",
370  // //20, t0_estim-100, t0_estim+100);
371  // 700, 0, 7000);
372  //}
373  //if(debug)
374  // cout<<" accepted + 400ns"<<endl;
375  //hT0SectorHisto->Fill((*lHisto).second->GetMean() + 400);
376  //}
377  if(debug)
378  cout<<endl;
379 
380  theT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetMean();
381  theSigmaT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetRMS();*/
382  }
383  /*if(!hT0SectorHisto){
384  cout<<"[DTT0CalibrationNew]: All the t0 per layer are still uncorrect: trying with greater number of events"<<endl;
385  eventsForLayerT0 = eventsForLayerT0*2;
386  return;
387  }
388  if(debug)
389  cout<<"[DTT0CalibrationNew] t0 reference for this sector "<<
390  hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())<<endl;*/
391  if(increaseEvtsForLayerT0){
392  cout<<"[DTT0CalibrationNew]: t0 per layer are still uncorrect: trying with greater number of events"<<endl;
394  return;
395  }
396  }
397 }
RunNumber_t run() const
Definition: EventID.h:39
EventNumber_t event() const
Definition: EventID.h:41
std::map< DTWireId, TH1I * > theHistoWireMap
DTChamberId chamberId() const
Return the corresponding ChamberId.
unsigned int eventsForLayerT0
unsigned int eventsForWireT0
unsigned int timeBoxWidth
DTSuperLayerId superlayerId() const
Return the corresponding SuperLayerId.
Definition: DTLayerId.h:59
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:20
std::map< DTLayerId, double > theTPPeakMap
std::map< DTChamberId, double > theSumT0ByChamber
std::map< DTChamberId, int > theCountT0ByChamber
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::map< DTWireId, int > nDigiPerWire
std::vector< DTWireId > wireIdWithHistos
std::map< DTWireId, double > mK
bin
set the eta bin as selection string.
std::map< DTWireId, int > nDigiPerWire_ref
std::map< DTLayerId, TH1I * > theHistoLayerMap
int get(int wheelId, int stationId, int sectorId, int slId, float &tTrig, float &tTrms, float &kFact, DTTimeUnits::type unit) const
get content
Definition: DTTtrig.cc:85
std::map< DTWireId, double > theAbsoluteT0PerWire
std::map< DTWireId, double > mK_ref
const T & get() const
Definition: EventSetup.h:56
std::vector< DTDigi >::const_iterator const_iterator
edm::EventID id() const
Definition: EventBase.h:60
int sector() const
Definition: DTChamberId.h:61
unsigned int rejectDigiFromPeak
std::pair< const_iterator, const_iterator > Range
std::string getHistoName(const DTWireId &wId) const
std::map< DTWireId, double > qK
std::string theCalibSector
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:45
edm::ESHandle< DTGeometry > dtGeom
void DTT0CalibrationNew::endJob ( void  )
virtual

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 399 of file DTT0CalibrationNew.cc.

References funct::abs(), DTT0::begin(), relativeConstraints::chamber, DTSuperLayerId::chamberId(), DTTimeUnits::counts, gather_cfg::cout, debug, DTT0::end(), DTT0::get(), DTT0::set(), mathSSE::sqrt(), cscNeutronWriter_cfi::t0, interactiveExample::theFile, tzero, and DTCalibDBUtils::writeToDB().

Referenced by o2olib.O2ORunMgr::executeJob().

399  {
400 
401  DTT0* t0s = new DTT0();
402  DTT0* t0sWRTChamber = new DTT0();
403 
404  if(debug)
405  cout << "[DTT0CalibrationNewPerLayer]Writing histos to file!" << endl;
406 
407  theFile->cd();
408  //hT0SectorHisto->Write();
409  for(map<DTWireId, TH1I*>::const_iterator wHisto = theHistoWireMap.begin();
410  wHisto != theHistoWireMap.end();
411  ++wHisto) {
412  (*wHisto).second->Write();
413  }
414  for(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin();
415  lHisto != theHistoLayerMap.end();
416  ++lHisto) {
417  (*lHisto).second->Write();
418  }
419 
420  if(debug)
421  cout << "[DTT0CalibrationNew] Compute and store t0 and sigma per wire" << endl;
422 
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]);
426 
427  for(map<DTWireId, double>::const_iterator wiret0 = theAbsoluteT0PerWire.begin();
428  wiret0 != theAbsoluteT0PerWire.end();
429  ++wiret0){
430  if(nDigiPerWire[(*wiret0).first]){
431  double t0 = (*wiret0).second/nDigiPerWire[(*wiret0).first];
432  DTChamberId chamberId = ((*wiret0).first).chamberId();
433  //theRelativeT0PerWire[(*wiret0).first] = t0 - hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin());
434  theRelativeT0PerWire[(*wiret0).first] = t0 - theRefT0ByChamber[chamberId];
435  cout<<"Wire "<<(*wiret0).first<<" has t0 "<<t0<<"(absolute) "<<theRelativeT0PerWire[(*wiret0).first]<<"(relative)";
436 
437  //theSigmaT0PerWire[(*wiret0).first] = sqrt((theSigmaT0PerWire[(*wiret0).first] / nDigiPerWire[(*wiret0).first]) - t0*t0);
438  theSigmaT0PerWire[(*wiret0).first] = sqrt(qK[(*wiret0).first]/nDigiPerWire[(*wiret0).first]);
439  cout<<" sigma "<<theSigmaT0PerWire[(*wiret0).first]<<endl;
440  }
441  else{
442  cout<<"[DTT0CalibrationNew] ERROR: no digis in wire "<<(*wiret0).first<<endl;
443  abort();
444  }
445  }
446 
448  // Get all the sls from the setup
449  const vector<const DTSuperLayer*>& superLayers = dtGeom->superLayers();
450  // Loop over all SLs
451  for(auto sl = superLayers.begin();
452  sl != superLayers.end(); sl++) {
453 
454 
455  //Compute mean for odd and even superlayers
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();
462  ++wiret0){
463  if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
464  if(debug)
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;
469  oddLayersDen++;
470  }
471  else{
472  evenLayersMean = evenLayersMean + (*wiret0).second;
473  evenLayersDen++;
474  }
475  }
476  }
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;
481 
482  //Compute sigma for odd and even superlayers
483  double oddLayersSigma=0;
484  double evenLayersSigma=0;
485  for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
486  wiret0 != theRelativeT0PerWire.end();
487  ++wiret0){
488  if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
489  if(((*wiret0).first.layerId().layer()) % 2){
490  oddLayersSigma = oddLayersSigma + ((*wiret0).second - oddLayersMean) * ((*wiret0).second - oddLayersMean);
491  }
492  else{
493  evenLayersSigma = evenLayersSigma + ((*wiret0).second - evenLayersMean) * ((*wiret0).second - evenLayersMean);
494  }
495  }
496  }
497  oddLayersSigma = oddLayersSigma/oddLayersDen;
498  evenLayersSigma = evenLayersSigma/evenLayersDen;
499  oddLayersSigma = sqrt(oddLayersSigma);
500  evenLayersSigma = sqrt(evenLayersSigma);
501 
502  if(debug && oddLayersMean)
503  cout<<"[DTT0CalibrationNew] Relative T0 sigma for odd layers "<<oddLayersSigma<<" even layers"<<evenLayersSigma<<endl;
504 
505  //Recompute the mean for odd and even superlayers discarding fluctations
506  double oddLayersFinalMean=0;
507  double evenLayersFinalMean=0;
508  for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
509  wiret0 != theRelativeT0PerWire.end();
510  ++wiret0){
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;
515  }
516  else{
517  if(abs((*wiret0).second - evenLayersMean) < (2*evenLayersSigma))
518  evenLayersFinalMean = evenLayersFinalMean + (*wiret0).second;
519  }
520  }
521  }
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;
526 
527  for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
528  wiret0 != theRelativeT0PerWire.end();
529  ++wiret0){
530  if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
531  double t0=-999;
532  if(((*wiret0).first.layerId().layer()) % 2)
533  t0 = (*wiret0).second + (evenLayersFinalMean - oddLayersFinalMean);
534  else
535  t0 = (*wiret0).second;
536 
537  cout<<"[DTT0CalibrationNew] Wire "<<(*wiret0).first<<" has t0 "<<(*wiret0).second<<" (relative, after even-odd layer corrections) "
538  <<" sigma "<<theSigmaT0PerWire[(*wiret0).first]<<endl;
539  //Store the results into DB
540  t0s->set((*wiret0).first, t0, theSigmaT0PerWire[(*wiret0).first],DTTimeUnits::counts);
541  }
542  }
543  }
544 
546  if(debug)
547  cout << "[DTT0CalibrationNew]Computing relative t0 wrt to chamber average" << endl;
548  //Compute the reference for each chamber
549  map<DTChamberId,double> sumT0ByChamber;
550  map<DTChamberId,int> countT0ByChamber;
551  for(DTT0::const_iterator tzero = t0s->begin();
552  tzero != t0s->end(); ++tzero) {
553 // @@@ NEW DTT0 FORMAT
554 // DTChamberId chamberId((*tzero).first.wheelId,
555 // (*tzero).first.stationId,
556 // (*tzero).first.sectorId);
557 // sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + (*tzero).second.t0mean;
558  int channelId = tzero->channelId;
559  if ( channelId == 0 ) continue;
560  DTWireId wireId(channelId);
561  DTChamberId chamberId(wireId.chamberId());
562  //sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + tzero->t0mean;
563 // @@@ better DTT0 usage
564  float t0mean_f;
565  float t0rms_f;
566  t0s->get(wireId,t0mean_f,t0rms_f,DTTimeUnits::counts);
567  sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + t0mean_f;
568 // @@@ NEW DTT0 END
569  countT0ByChamber[chamberId]++;
570  }
571 
572  //Change reference for each wire and store the new t0s in the new map
573  for(DTT0::const_iterator tzero = t0s->begin();
574  tzero != t0s->end(); ++tzero) {
575 // @@@ NEW DTT0 FORMAT
576 // DTChamberId chamberId((*tzero).first.wheelId,
577 // (*tzero).first.stationId,
578 // (*tzero).first.sectorId);
579 // double t0mean = ((*tzero).second.t0mean) - (sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
580 // double t0rms = (*tzero).second.t0rms;
581 // DTWireId wireId((*tzero).first.wheelId,
582 // (*tzero).first.stationId,
583 // (*tzero).first.sectorId,
584 // (*tzero).first.slId,
585 // (*tzero).first.layerId,
586 // (*tzero).first.cellId);
587  int channelId = tzero->channelId;
588  if ( channelId == 0 ) continue;
589  DTWireId wireId( channelId );
590  DTChamberId chamberId(wireId.chamberId());
591  //double t0mean = (tzero->t0mean) - (sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
592  //double t0rms = tzero->t0rms;
593 // @@@ better DTT0 usage
594  float t0mean_f;
595  float t0rms_f;
596  t0s->get(wireId,t0mean_f,t0rms_f,DTTimeUnits::counts);
597  double t0mean = t0mean_f - (sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
598  double t0rms = t0rms_f;
599 // @@@ NEW DTT0 END
600  t0sWRTChamber->set(wireId,
601  t0mean,
602  t0rms,
604  if(debug){
605  //cout<<"Chamber "<<chamberId<<" has reference "<<(sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
606 // cout<<"Changing t0 of wire "<<wireId<<" from "<<(*tzero).second.t0mean<<" to "<<t0mean<<endl;
607  cout<<"Changing t0 of wire "<<wireId<<" from "<<tzero->t0mean<<" to "<<t0mean<<endl;
608  }
609  }
610 
612  if(debug)
613  cout << "[DTT0CalibrationNew]Writing values in DB!" << endl;
614  // FIXME: to be read from cfg?
615  string t0Record = "DTT0Rcd";
616  // Write the t0 map to DB
617  DTCalibDBUtils::writeToDB(t0Record, t0sWRTChamber);
618 }
const_iterator begin() const
Definition: DTT0.cc:204
int set(int wheelId, int stationId, int sectorId, int slId, int layerId, int cellId, float t0mean, float t0rms, DTTimeUnits::type unit)
Definition: DTT0.cc:140
std::map< DTChamberId, double > theRefT0ByChamber
std::map< DTWireId, TH1I * > theHistoWireMap
int get(int wheelId, int stationId, int sectorId, int slId, int layerId, int cellId, float &t0mean, float &t0rms, DTTimeUnits::type unit) const
Definition: DTT0.cc:67
std::vector< DTT0Data >::const_iterator const_iterator
Access methods to data.
Definition: DTT0.h:140
Definition: DTT0.h:53
std::map< DTChamberId, double > theSumT0ByChamber
std::map< DTWireId, double > theSigmaT0PerWire
T sqrt(T t)
Definition: SSEVec.h:18
std::map< DTChamberId, int > theCountT0ByChamber
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::map< DTWireId, int > nDigiPerWire
std::map< DTWireId, double > theRelativeT0PerWire
const_iterator end() const
Definition: DTT0.cc:209
std::map< DTLayerId, TH1I * > theHistoLayerMap
std::map< DTWireId, double > theAbsoluteT0PerWire
static const double tzero[3]
std::map< DTWireId, double > qK
edm::ESHandle< DTGeometry > dtGeom
static void writeToDB(std::string record, T *payload)
const std::vector< const DTSuperLayer * > & superLayers() const
Return a vector of all SuperLayer.
Definition: DTGeometry.cc:90
string DTT0CalibrationNew::getHistoName ( const DTWireId wId) const
private

Definition at line 620 of file DTT0CalibrationNew.cc.

References DTLayerId::layer(), DTChamberId::sector(), DTChamberId::station(), DTSuperLayerId::superlayer(), DTChamberId::wheel(), and DTWireId::wire().

620  {
621  string histoName;
622  stringstream theStream;
623  theStream << "Ch_" << wId.wheel() << "_" << wId.station() << "_" << wId.sector()
624  << "_SL" << wId.superlayer() << "_L" << wId.layer() << "_W"<< wId.wire() <<"_hT0Histo";
625  theStream >> histoName;
626  return histoName;
627 }
int layer() const
Return the layer number.
Definition: DTLayerId.h:53
int wire() const
Return the wire number.
Definition: DTWireId.h:56
int superlayer() const
Return the superlayer number (deprecated method name)
int sector() const
Definition: DTChamberId.h:61
int station() const
Return the station number.
Definition: DTChamberId.h:51
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:45
string DTT0CalibrationNew::getHistoName ( const DTLayerId lId) const
private

Definition at line 629 of file DTT0CalibrationNew.cc.

References DTLayerId::layer(), DTChamberId::sector(), DTChamberId::station(), DTSuperLayerId::superlayer(), and DTChamberId::wheel().

629  {
630  string histoName;
631  stringstream theStream;
632  theStream << "Ch_" << lId.wheel() << "_" << lId.station() << "_" << lId.sector()
633  << "_SL" << lId.superlayer() << "_L" << lId.layer() <<"_hT0Histo";
634  theStream >> histoName;
635  return histoName;
636 }
int layer() const
Return the layer number.
Definition: DTLayerId.h:53
int superlayer() const
Return the superlayer number (deprecated method name)
int sector() const
Definition: DTChamberId.h:61
int station() const
Return the station number.
Definition: DTChamberId.h:51
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:45

Member Data Documentation

std::vector<std::string> DTT0CalibrationNew::cellsWithHistos
private

Definition at line 102 of file DTT0CalibrationNew.h.

std::string DTT0CalibrationNew::dbLabel
private

Definition at line 59 of file DTT0CalibrationNew.h.

bool DTT0CalibrationNew::debug
private
std::string DTT0CalibrationNew::digiLabel
private

Definition at line 57 of file DTT0CalibrationNew.h.

edm::ESHandle<DTGeometry> DTT0CalibrationNew::dtGeom
private

Definition at line 125 of file DTT0CalibrationNew.h.

unsigned int DTT0CalibrationNew::eventsForLayerT0
private

Definition at line 69 of file DTT0CalibrationNew.h.

unsigned int DTT0CalibrationNew::eventsForWireT0
private

Definition at line 71 of file DTT0CalibrationNew.h.

TH1D* DTT0CalibrationNew::hT0SectorHisto
private

Definition at line 96 of file DTT0CalibrationNew.h.

std::map<DTWireId,double> DTT0CalibrationNew::mK
private

Definition at line 110 of file DTT0CalibrationNew.h.

std::map<DTWireId,double> DTT0CalibrationNew::mK_ref
private

Definition at line 111 of file DTT0CalibrationNew.h.

std::map<DTWireId,int> DTT0CalibrationNew::nDigiPerWire
private

Definition at line 108 of file DTT0CalibrationNew.h.

std::map<DTWireId,int> DTT0CalibrationNew::nDigiPerWire_ref
private

Definition at line 109 of file DTT0CalibrationNew.h.

unsigned int DTT0CalibrationNew::nevents
private

Definition at line 67 of file DTT0CalibrationNew.h.

std::map<DTWireId,double> DTT0CalibrationNew::qK
private

Definition at line 112 of file DTT0CalibrationNew.h.

unsigned int DTT0CalibrationNew::rejectDigiFromPeak
private

Definition at line 83 of file DTT0CalibrationNew.h.

unsigned int DTT0CalibrationNew::retryForLayerT0
private

Definition at line 85 of file DTT0CalibrationNew.h.

int DTT0CalibrationNew::selSector
private

Definition at line 91 of file DTT0CalibrationNew.h.

int DTT0CalibrationNew::selWheel
private

Definition at line 89 of file DTT0CalibrationNew.h.

TSpectrum* DTT0CalibrationNew::spectrum
private

Definition at line 98 of file DTT0CalibrationNew.h.

std::map<DTWireId,double> DTT0CalibrationNew::theAbsoluteT0PerWire
private

Definition at line 105 of file DTT0CalibrationNew.h.

std::string DTT0CalibrationNew::theCalibSector
private

Definition at line 90 of file DTT0CalibrationNew.h.

std::string DTT0CalibrationNew::theCalibWheel
private

Definition at line 88 of file DTT0CalibrationNew.h.

std::map<DTChamberId,int> DTT0CalibrationNew::theCountT0ByChamber
private

Definition at line 121 of file DTT0CalibrationNew.h.

TFile* DTT0CalibrationNew::theFile
private

Definition at line 62 of file DTT0CalibrationNew.h.

std::map<DTLayerId, TH1I*> DTT0CalibrationNew::theHistoLayerMap
private

Definition at line 94 of file DTT0CalibrationNew.h.

std::map<DTWireId,TH1I*> DTT0CalibrationNew::theHistoWireMap
private

Definition at line 114 of file DTT0CalibrationNew.h.

TFile* DTT0CalibrationNew::theOutputFile
private

Definition at line 64 of file DTT0CalibrationNew.h.

std::map<DTChamberId,double> DTT0CalibrationNew::theRefT0ByChamber
private

Definition at line 122 of file DTT0CalibrationNew.h.

std::map<DTWireId,double> DTT0CalibrationNew::theRelativeT0PerWire
private

Definition at line 106 of file DTT0CalibrationNew.h.

std::map<std::string,double> DTT0CalibrationNew::theSigmaT0LayerMap
private

Definition at line 117 of file DTT0CalibrationNew.h.

std::map<DTWireId,double> DTT0CalibrationNew::theSigmaT0PerWire
private

Definition at line 107 of file DTT0CalibrationNew.h.

std::map<DTChamberId,double> DTT0CalibrationNew::theSumT0ByChamber
private

Definition at line 120 of file DTT0CalibrationNew.h.

std::map<std::string,double> DTT0CalibrationNew::theT0LayerMap
private

Definition at line 116 of file DTT0CalibrationNew.h.

std::map<DTLayerId,double> DTT0CalibrationNew::theTPPeakMap
private

Definition at line 118 of file DTT0CalibrationNew.h.

unsigned int DTT0CalibrationNew::timeBoxWidth
private

Definition at line 80 of file DTT0CalibrationNew.h.

double DTT0CalibrationNew::tpPeakWidth
private

Definition at line 74 of file DTT0CalibrationNew.h.

double DTT0CalibrationNew::tpPeakWidthPerLayer
private

Definition at line 77 of file DTT0CalibrationNew.h.

std::vector<DTWireId> DTT0CalibrationNew::wireIdWithHistos
private

Definition at line 101 of file DTT0CalibrationNew.h.