#include <DTT0CalibrationNew.h>

Inheritance diagram for DTT0CalibrationNew:

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< ConsumesInfo >	consumesInfo () const

	EDConsumerBase ()

ProductHolderIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const

std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesDependentUpon (std::string const &iProcessName, std::string const &iModuleLabel, bool iPrint, std::vector< char const * > &oModuleLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

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< 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

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, dtTTrigAnalyzer_cfg::dbLabel, debug, dtTPAnalyzer_cfg::digiLabel, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), cmsHarvester::nevents, dtT0Analyzer_cfg::rootFileName, relativeConstraints::station, and interactiveExample::theFile.

                                                                   {
   // Get the debug parameter for verbose output
   debug = pset.getUntrackedParameter<bool>("debug");
   if(debug) 
     cout << "[DTT0CalibrationNew]Constructor called!" << endl;
 
   // Get the label to retrieve digis from the event
   digiLabel = pset.getUntrackedParameter<string>("digiLabel");
 
   dbLabel  = pset.getUntrackedParameter<string>("dbLabel", "");
 
   // The root file which contain the histos per layer
   string rootFileName = pset.getUntrackedParameter<string>("rootFileName","DTT0PerLayer.root");
   theFile = new TFile(rootFileName.c_str(), "RECREATE");
  
   theCalibWheel =  pset.getUntrackedParameter<string>("calibWheel", "All"); //FIXME amke a vector of integer instead of a string
   if(theCalibWheel != "All") {
     stringstream linestr;
     int selWheel;
     linestr << theCalibWheel;
     linestr >> selWheel;
     cout << "[DTT0CalibrationNewPerLayer] chosen wheel " << selWheel << endl;
   }
 
   // Sector/s to calibrate
   theCalibSector =  pset.getUntrackedParameter<string>("calibSector", "All"); //FIXME amke a vector of integer instead of a string
   if(theCalibSector != "All") {
     stringstream linestr;
     int selSector;
     linestr << theCalibSector;
     linestr >> selSector;
     cout << "[DTT0CalibrationNewPerLayer] chosen sector " << selSector << endl;
   }
 
   vector<string> defaultCell;
   defaultCell.push_back("None");
   cellsWithHistos = pset.getUntrackedParameter<vector<string> >("cellsWithHisto", defaultCell);
   for(vector<string>::const_iterator cell = cellsWithHistos.begin(); cell != cellsWithHistos.end(); ++cell){
     if((*cell) != "None"){
       stringstream linestr;
       int wheel,sector,station,sl,layer,wire;
       linestr << (*cell);
       linestr >> wheel >> sector >> station >> sl >> layer >> wire;
       wireIdWithHistos.push_back(DTWireId(wheel,station,sector,sl,layer,wire));
     }
   }
 
   hT0SectorHisto=0;
 
   nevents=0;
   eventsForLayerT0 = pset.getParameter<unsigned int>("eventsForLayerT0");
   eventsForWireT0 = pset.getParameter<unsigned int>("eventsForWireT0");
   tpPeakWidth = pset.getParameter<double>("tpPeakWidth");
   tpPeakWidthPerLayer = pset.getParameter<double>("tpPeakWidthPerLayer");
   timeBoxWidth = pset.getParameter<unsigned int>("timeBoxWidth"); 
   rejectDigiFromPeak = pset.getParameter<unsigned int>("rejectDigiFromPeak"); 
 
   spectrum = new TSpectrum(5);
   retryForLayerT0 = 0;  
 }

DTT0CalibrationNew::~DTT0CalibrationNew ( )

virtual

Destructor.

Definition at line 96 of file DTT0CalibrationNew.cc.

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

                                        {
   if(debug) 
     cout << "[DTT0CalibrationNew]Destructor called!" << endl;
 
   delete spectrum;
   theFile->Close();
 }

Member Function Documentation

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(), newFWLiteAna::bin, DTSuperLayerId::chamberId(), DTTimeUnits::counts, gather_cfg::cout, dtTTrigAnalyzer_cfg::dbLabel, debug, dtTPAnalyzer_cfg::digiLabel, edm::EventID::event(), spr::find(), edm::EventSetup::get(), mergeVDriftHistosByStation::getHistoName(), edm::EventBase::id(), timingPdfMaker::mean, cmsHarvester::nevents, edm::EventID::run(), DTChamberId::sector(), DTLayerId::superlayerId(), interactiveExample::theFile, and DTChamberId::wheel().

                                                                                         {
   if(debug || event.id().event() % 500==0)
     cout << "--- [DTT0CalibrationNew] Analysing Event: #Run: " << event.id().run()
      << " #Event: " << event.id().event() << endl;
   nevents++;
 
   // Get the digis from the event
   Handle<DTDigiCollection> digis; 
   event.getByLabel(digiLabel, digis);
 
   // Get the DT Geometry
   eventSetup.get<MuonGeometryRecord>().get(dtGeom);
 
   // Get ttrig DB
   edm::ESHandle<DTTtrig> tTrigMap;
   eventSetup.get<DTTtrigRcd>().get(dbLabel,tTrigMap);
   
   // Iterate through all digi collections ordered by LayerId   
   DTDigiCollection::DigiRangeIterator dtLayerIt;
   for (dtLayerIt = digis->begin();
        dtLayerIt != digis->end();
        ++dtLayerIt){
     // Get the iterators over the digis associated with this LayerId
     const DTDigiCollection::Range& digiRange = (*dtLayerIt).second;
   
     // Get the layerId
     const DTLayerId layerId = (*dtLayerIt).first; //FIXME: check to be in the right sector
     const DTChamberId chamberId = layerId.superlayerId().chamberId();
 
     if((theCalibWheel != "All") && (layerId.superlayerId().chamberId().wheel() != selWheel))
       continue;
     if((theCalibSector != "All") && (layerId.superlayerId().chamberId().sector() != selSector))
       continue;
  
     //if(debug) {
     //  cout << "Layer " << layerId<<" with "<<distance(digiRange.first, digiRange.second)<<" digi"<<endl;
     //}
 
     float tTrig,tTrigRMS, kFactor;
     tTrigMap->get(layerId.superlayerId(), tTrig, tTrigRMS, kFactor, DTTimeUnits::counts );
     if(debug&&(nevents <= 1)){
     cout << "  Superlayer: " << layerId.superlayerId() << endl 
          << "            tTrig,tTrigRMS= " << tTrig << ", " << tTrigRMS << endl;
     }   
 
     // Loop over all digis in the given layer
     for (DTDigiCollection::const_iterator digi = digiRange.first;
      digi != digiRange.second;
      ++digi) {
       double t0 = (*digi).countsTDC();
 
       //Use first bunch of events to fill t0 per layer
       if(nevents < eventsForLayerT0){
     //Get the per-layer histo from the map
     TH1I *hT0LayerHisto = theHistoLayerMap[layerId];
     //If it doesn't exist, book it
     if(hT0LayerHisto == 0){
       theFile->cd();
       float hT0Min = tTrig - 2*tTrigRMS;
       float hT0Max = hT0Min + timeBoxWidth;
       hT0LayerHisto = new TH1I(getHistoName(layerId).c_str(),
                    "T0 from pulses by layer (TDC counts, 1 TDC count = 0.781 ns)",
                    timeBoxWidth,hT0Min,hT0Max);
       if(debug)
         cout << "  New T0 per Layer Histo: " << hT0LayerHisto->GetName() << endl;
       theHistoLayerMap[layerId] = hT0LayerHisto;
     }
     
     //Fill the histos
     theFile->cd();
     if(hT0LayerHisto != 0) {
       //  if(debug)
       // cout<<"Filling histo "<<hT0LayerHisto->GetName()<<" with digi "<<t0<<" TDC counts"<<endl;
       hT0LayerHisto->Fill(t0);
     }
       }
 
       //Use all the remaining events to compute t0 per wire
       if(nevents>eventsForLayerT0){
     // Get the wireId
     const DTWireId wireId(layerId, (*digi).wire());
     if(debug) {
       cout << "   Wire: " << wireId << endl
            << "       time (TDC counts): " << (*digi).countsTDC()<< endl;
     }   
 
     //Fill the histos per wire for the chosen cells
     vector<DTWireId>::iterator it = find(wireIdWithHistos.begin(),wireIdWithHistos.end(),wireId);
     if (it!=wireIdWithHistos.end()){
       //Get the per-wire histo from the map
       TH1I *hT0WireHisto = theHistoWireMap[wireId]; 
       //If it doesn't exist, book it
       if(hT0WireHisto == 0){
         theFile->cd(); 
         hT0WireHisto = new TH1I(getHistoName(wireId).c_str(),"T0 from pulses by wire (TDC counts, 1 TDC count = 0.781 ns)",7000,0,7000);
         //hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())-100,
         //hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())+100);
         if(debug)
           cout << "  New T0 per wire Histo: " << hT0WireHisto->GetName() << endl;
         theHistoWireMap[wireId] = hT0WireHisto;
       }
       //Fill the histos
       theFile->cd();
       if(hT0WireHisto != 0) {
         //if(debug)
         // cout<<"Filling histo "<<hT0WireHisto->GetName()<<" with digi "<<t0<<" TDC counts"<<endl;
         hT0WireHisto->Fill(t0);
       }
     }
 
     //Check the tzero has reasonable value
     //float hT0Min = tTrig - 2*tTrigRMS;
     //float hT0Max = hT0Min + timeBoxWidth;
     /*if(abs(hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin()) - t0) > rejectDigiFromPeak){
       if(debug)
         cout<<"digi skipped because t0 per sector "<<hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())<<endl;
       continue;
     }*/
     /*if((t0 < hT0Min)||(t0 > hT0Max)){
           if(debug)
             cout<<"digi skipped because t0 outside of interval (" << hT0Min << "," << hT0Max << ")" <<endl;
           continue;
         }*/
     //Select per layer
     if(fabs(theTPPeakMap[layerId] - t0) > rejectDigiFromPeak){
       if(debug)
             cout<<"digi skipped because t0 too far from peak " << theTPPeakMap[layerId] << endl;
           continue; 
     }
 
     //Find to ref. per chamber
     theSumT0ByChamber[chamberId] = theSumT0ByChamber[chamberId] + t0;
     theCountT0ByChamber[chamberId]++;
 
     //Use second bunch of events to compute a t0 reference per wire
     if(nevents< (eventsForLayerT0 + eventsForWireT0)){
       if(!nDigiPerWire_ref[wireId]){
         mK_ref[wireId] = 0;
       }
       nDigiPerWire_ref[wireId] = nDigiPerWire_ref[wireId] + 1;
       mK_ref[wireId] = mK_ref[wireId] + (t0-mK_ref[wireId])/nDigiPerWire_ref[wireId];
     }
     //Use last all the remaining events to compute the mean and sigma t0 per wire
     else if(nevents>(eventsForLayerT0 + eventsForWireT0)){
       if(abs(t0-mK_ref[wireId]) > tpPeakWidth)
         continue;
       if(!nDigiPerWire[wireId]){
         theAbsoluteT0PerWire[wireId] = 0;
         qK[wireId] = 0;
         mK[wireId] = 0;
       }
       nDigiPerWire[wireId] = nDigiPerWire[wireId] + 1;
       theAbsoluteT0PerWire[wireId] = theAbsoluteT0PerWire[wireId] + t0;
       //theSigmaT0PerWire[wireId] = theSigmaT0PerWire[wireId] + (t0*t0);
       qK[wireId] = qK[wireId] + ((nDigiPerWire[wireId]-1)*(t0-mK[wireId])*(t0-mK[wireId])/nDigiPerWire[wireId]);
       mK[wireId] = mK[wireId] + (t0-mK[wireId])/nDigiPerWire[wireId];
     }
       }//end if(nevents>1000)
     }//end loop on digi
   }//end loop on layer
 
   //Use the t0 per layer histos to have an indication about the t0 position 
   if(nevents == eventsForLayerT0){
     bool increaseEvtsForLayerT0 = false;    
     for(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin();
     lHisto != theHistoLayerMap.end();
     ++lHisto){
       if(debug)
     cout<<"Reading histogram "<<(*lHisto).second->GetName()<<" with mean "<<(*lHisto).second->GetMean()<<" and RMS "<<(*lHisto).second->GetRMS() << endl;
 
       //Find peaks
       //int npeaks = spectrum->Search((*lHisto).second,0.5,"goff");
       //int npeaks = spectrum->Search((*lHisto).second,(tpPeakWidthPerLayer/2.),"goff",0.3);
       int npeaks = spectrum->Search((*lHisto).second,(tpPeakWidthPerLayer/2.),"",0.3);
 
       double *peaks = spectrum->GetPositionX();
       //Put in a std::vector<float>
       vector<float> peakMeans(peaks,peaks + npeaks);
       //Sort the peaks in ascending order
       sort(peakMeans.begin(),peakMeans.end());
                 
       //Find best peak -- preliminary criteria: find peak closest to center of time box 
       float tTrig,tTrigRMS, kFactor;
       tTrigMap->get((*lHisto).first.superlayerId(), tTrig, tTrigRMS, kFactor, DTTimeUnits::counts );
 
       float timeBoxCenter = (2*tTrig + (float)timeBoxWidth)/2.; 
       float hMin = (*lHisto).second->GetXaxis()->GetXmin();
       float hMax = (*lHisto).second->GetXaxis()->GetXmax();     
       vector<float>::const_iterator tpPeak = peakMeans.end();
       for(vector<float>::const_iterator it = peakMeans.begin(); it != peakMeans.end(); ++it){
     float mean = *it;
 
     int bin = (*lHisto).second->GetXaxis()->FindBin(mean);
     float yp = (*lHisto).second->GetBinContent(bin);
     if(debug) cout << "Peak : (" << mean << "," << yp << ")" << endl; 
 
     //Find RMS
     float previous_peak = (it == peakMeans.begin())?hMin:*(it - 1);
         float next_peak = (it == (peakMeans.end()-1))?hMax:*(it + 1);
 
     float rangemin = mean - (mean - previous_peak)/8.;
         float rangemax = mean + (next_peak - mean)/8.;
     int binmin = (*lHisto).second->GetXaxis()->FindBin(rangemin);
     int binmax = (*lHisto).second->GetXaxis()->FindBin(rangemax);
     (*lHisto).second->GetXaxis()->SetRange(binmin,binmax);
     //RMS estimate
     float rms_seed = (*lHisto).second->GetRMS();
 
     /*rangemin = mean - 2*rms_seed;
     rangemax = mean + 2*rms_seed;
     if(debug) cout << "Seed for RMS, Fit min, Fit max: " << rms_seed << ", " << rangemin << ", " << rangemax << endl;
     //Fit to gaussian
     string funcname("fitFcn_");
     funcname += (*lHisto).second->GetName();
     if(debug) cout << "Fitting function " << funcname << endl; 
     TF1* func = new TF1(funcname.c_str(),"gaus",rangemin,rangemax);
     func->SetParameters(yp,mean,rms_seed);
     (*lHisto).second->Fit(func,"Q","",rangemin,rangemax);
 
     float fitconst = func->GetParameter(0);
     float fitmean = func->GetParameter(1);
     float fitrms = func->GetParameter(2);
     float chisquare = func->GetChisquare()/func->GetNDF();
     if(debug) cout << "Gaussian fit constant,mean,RMS,chi2= " << fitconst << ", " << fitmean << ", " << fitrms << ", " << chisquare << endl;*/
 
     //Reject peaks with RMS larger than specified
     //if(fitrms > tpPeakWidth) continue;
     if(rms_seed > tpPeakWidthPerLayer) continue;
 
     if(fabs(mean - timeBoxCenter) < fabs(*tpPeak - timeBoxCenter)) tpPeak = it;
       } 
       //Didn't find peak    
       /*if(tpPeak == peakMeans.end()){
     if(retryForLayerT0 < 2){
       increaseEvtsForLayerT0 = true;
       retryForLayerT0++;
       break;
     } 
       }*/
 
       float selPeak = (tpPeak != peakMeans.end())?*tpPeak:(*lHisto).second->GetBinCenter((*lHisto).second->GetMaximumBin());        
       if(debug) cout << "Peak selected at " << selPeak << endl;
     
       theTPPeakMap[(*lHisto).first] = selPeak;
         
       //Take the mean as a first t0 estimation
       /*if((*lHisto).second->GetRMS() < tpPeakWidth){
     if(hT0SectorHisto == 0){
       hT0SectorHisto = new TH1D("hT0AllLayerOfSector","T0 from pulses per layer in sector", 
                     //20, (*lHisto).second->GetMean()-100, (*lHisto).second->GetMean()+100);
                     700, 0, 7000);
                     //300,3300,3600);   
     }
     if(debug)
       cout<<" accepted"<<endl;
     //TH1I* aux_histo = (*lHisto).second;
     //aux_histo->GetXaxis()->SetRangeUser(3300,3600);
     hT0SectorHisto->Fill((*lHisto).second->GetMean());
     //hT0SectorHisto->Fill(aux_histo->GetMean());
       }
       //Take the mean of noise + 400ns as a first t0 estimation
       //if((*lHisto).second->GetRMS()>10.0 && ((*lHisto).second->GetRMS()<15.0)){
       //double t0_estim = (*lHisto).second->GetMean() + 400;
       //if(hT0SectorHisto == 0){
       //  hT0SectorHisto = new TH1D("hT0AllLayerOfSector","T0 from pulses per layer in sector", 
       //                //20, t0_estim-100, t0_estim+100);
       //                700, 0, 7000);
       //}
       //if(debug)
       //  cout<<" accepted + 400ns"<<endl;
       //hT0SectorHisto->Fill((*lHisto).second->GetMean() + 400);
       //}
       if(debug)
     cout<<endl;
 
       theT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetMean();
       theSigmaT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetRMS();*/
     }
     /*if(!hT0SectorHisto){
       cout<<"[DTT0CalibrationNew]: All the t0 per layer are still uncorrect: trying with greater number of events"<<endl;
       eventsForLayerT0 = eventsForLayerT0*2;
       return;
     }
     if(debug)
       cout<<"[DTT0CalibrationNew] t0 reference for this sector "<<
     hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())<<endl;*/
     if(increaseEvtsForLayerT0){
     cout<<"[DTT0CalibrationNew]: t0 per layer are still uncorrect: trying with greater number of events"<<endl;
     eventsForLayerT0 = eventsForLayerT0*2;
     return;
     }       
   } 
 }

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(), DTSuperLayerId::chamberId(), DTTimeUnits::counts, gather_cfg::cout, debug, DTT0::end(), DTT0::get(), DTT0::set(), mathSSE::sqrt(), interactiveExample::theFile, tzero, and DTCalibDBUtils::writeToDB().

Referenced by o2o.O2ORunMgr::executeJob().

                                 {
 
   DTT0* t0s = new DTT0();
   DTT0* t0sWRTChamber = new DTT0();
 
   if(debug) 
     cout << "[DTT0CalibrationNewPerLayer]Writing histos to file!" << endl;
 
   theFile->cd();
   //hT0SectorHisto->Write();
   for(map<DTWireId, TH1I*>::const_iterator wHisto = theHistoWireMap.begin();
       wHisto != theHistoWireMap.end();
       ++wHisto) {
     (*wHisto).second->Write(); 
   }
   for(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin();
       lHisto != theHistoLayerMap.end();
       ++lHisto) {
     (*lHisto).second->Write(); 
   }  
 
   if(debug) 
     cout << "[DTT0CalibrationNew] Compute and store t0 and sigma per wire" << endl;
 
   for(map<DTChamberId,double>::const_iterator chamber = theSumT0ByChamber.begin();
       chamber != theSumT0ByChamber.end();
       ++chamber) theRefT0ByChamber[(*chamber).first] = theSumT0ByChamber[(*chamber).first]/((double)theCountT0ByChamber[(*chamber).first]);
   
   for(map<DTWireId, double>::const_iterator wiret0 = theAbsoluteT0PerWire.begin();
       wiret0 != theAbsoluteT0PerWire.end();
       ++wiret0){
     if(nDigiPerWire[(*wiret0).first]){
       double t0 = (*wiret0).second/nDigiPerWire[(*wiret0).first];
       DTChamberId chamberId = ((*wiret0).first).chamberId();
       //theRelativeT0PerWire[(*wiret0).first] = t0 - hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin());
       theRelativeT0PerWire[(*wiret0).first] = t0 - theRefT0ByChamber[chamberId];    
       cout<<"Wire "<<(*wiret0).first<<" has    t0 "<<t0<<"(absolute) "<<theRelativeT0PerWire[(*wiret0).first]<<"(relative)";
 
       //theSigmaT0PerWire[(*wiret0).first] = sqrt((theSigmaT0PerWire[(*wiret0).first] / nDigiPerWire[(*wiret0).first]) - t0*t0);
       theSigmaT0PerWire[(*wiret0).first] = sqrt(qK[(*wiret0).first]/nDigiPerWire[(*wiret0).first]);
       cout<<"    sigma "<<theSigmaT0PerWire[(*wiret0).first]<<endl;
     }
     else{
       cout<<"[DTT0CalibrationNew] ERROR: no digis in wire "<<(*wiret0).first<<endl;
       abort();
     }
   }
 
   // Get all the sls from the setup
   const vector<const DTSuperLayer*>& superLayers = dtGeom->superLayers();     
   // Loop over all SLs
   for(auto  sl = superLayers.begin();
       sl != superLayers.end(); sl++) {
 
   
     //Compute mean for odd and even superlayers
     double oddLayersMean=0;
     double evenLayersMean=0; 
     double oddLayersDen=0;
     double evenLayersDen=0;
     for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
     wiret0 != theRelativeT0PerWire.end();
     ++wiret0){
       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
     if(debug)
       cout<<"[DTT0CalibrationNew] Superlayer "<<(*sl)->id()
           <<"layer " <<(*wiret0).first.layerId().layer()<<" with "<<(*wiret0).second<<endl;
     if(((*wiret0).first.layerId().layer()) % 2){
       oddLayersMean = oddLayersMean + (*wiret0).second;
       oddLayersDen++;
     }
     else{
       evenLayersMean = evenLayersMean + (*wiret0).second;
       evenLayersDen++;
     }
       }
     }
     oddLayersMean = oddLayersMean/oddLayersDen;
     evenLayersMean = evenLayersMean/evenLayersDen;
     if(debug && oddLayersMean)
       cout<<"[DTT0CalibrationNew] Relative T0 mean for  odd layers "<<oddLayersMean<<"  even layers"<<evenLayersMean<<endl;
 
     //Compute sigma for odd and even superlayers
     double oddLayersSigma=0;
     double evenLayersSigma=0;
     for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
     wiret0 != theRelativeT0PerWire.end();
     ++wiret0){
       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
     if(((*wiret0).first.layerId().layer()) % 2){
       oddLayersSigma = oddLayersSigma + ((*wiret0).second - oddLayersMean) * ((*wiret0).second - oddLayersMean);
     }
     else{
       evenLayersSigma = evenLayersSigma + ((*wiret0).second - evenLayersMean) * ((*wiret0).second - evenLayersMean);
     }
       }
     }
     oddLayersSigma = oddLayersSigma/oddLayersDen;
     evenLayersSigma = evenLayersSigma/evenLayersDen;
     oddLayersSigma = sqrt(oddLayersSigma);
     evenLayersSigma = sqrt(evenLayersSigma);
 
     if(debug && oddLayersMean)
       cout<<"[DTT0CalibrationNew] Relative T0 sigma for  odd layers "<<oddLayersSigma<<"  even layers"<<evenLayersSigma<<endl;
 
     //Recompute the mean for odd and even superlayers discarding fluctations
     double oddLayersFinalMean=0; 
     double evenLayersFinalMean=0;
     for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
     wiret0 != theRelativeT0PerWire.end();
     ++wiret0){
       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
     if(((*wiret0).first.layerId().layer()) % 2){
       if(abs((*wiret0).second - oddLayersMean) < (2*oddLayersSigma))
         oddLayersFinalMean = oddLayersFinalMean + (*wiret0).second;
     }
     else{
       if(abs((*wiret0).second - evenLayersMean) < (2*evenLayersSigma))
         evenLayersFinalMean = evenLayersFinalMean + (*wiret0).second;
     }
       }
     }
     oddLayersFinalMean = oddLayersFinalMean/oddLayersDen;
     evenLayersFinalMean = evenLayersFinalMean/evenLayersDen;
     if(debug && oddLayersMean)
       cout<<"[DTT0CalibrationNew] Final relative T0 mean for  odd layers "<<oddLayersFinalMean<<"  even layers"<<evenLayersFinalMean<<endl;
 
     for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
     wiret0 != theRelativeT0PerWire.end();
     ++wiret0){
       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
     double t0=-999;
     if(((*wiret0).first.layerId().layer()) % 2)
       t0 = (*wiret0).second + (evenLayersFinalMean - oddLayersFinalMean);
     else
       t0 = (*wiret0).second;
 
     cout<<"[DTT0CalibrationNew] Wire "<<(*wiret0).first<<" has    t0 "<<(*wiret0).second<<" (relative, after even-odd layer corrections)  "
         <<"    sigma "<<theSigmaT0PerWire[(*wiret0).first]<<endl;
     //Store the results into DB
     t0s->set((*wiret0).first, t0, theSigmaT0PerWire[(*wiret0).first],DTTimeUnits::counts); 
       }
     }
   }
   
   if(debug) 
     cout << "[DTT0CalibrationNew]Computing relative t0 wrt to chamber average" << endl;
   //Compute the reference for each chamber
   map<DTChamberId,double> sumT0ByChamber;
   map<DTChamberId,int> countT0ByChamber;
   for(DTT0::const_iterator tzero = t0s->begin();
       tzero != t0s->end(); ++tzero) {
 // @@@ NEW DTT0 FORMAT
 //    DTChamberId chamberId((*tzero).first.wheelId,
 //            (*tzero).first.stationId,
 //            (*tzero).first.sectorId);
 //    sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + (*tzero).second.t0mean;
     int channelId = tzero->channelId;
     if ( channelId == 0 ) continue;
     DTWireId wireId(channelId);
     DTChamberId chamberId(wireId.chamberId());
     //sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + tzero->t0mean;
 // @@@ better DTT0 usage
     float t0mean_f;
     float t0rms_f;
     t0s->get(wireId,t0mean_f,t0rms_f,DTTimeUnits::counts);
     sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + t0mean_f;
 // @@@ NEW DTT0 END
     countT0ByChamber[chamberId]++;
   }
 
   //Change reference for each wire and store the new t0s in the new map
   for(DTT0::const_iterator tzero = t0s->begin();
       tzero != t0s->end(); ++tzero) {
 // @@@ NEW DTT0 FORMAT
 //    DTChamberId chamberId((*tzero).first.wheelId,
 //            (*tzero).first.stationId,
 //            (*tzero).first.sectorId);
 //    double t0mean = ((*tzero).second.t0mean) - (sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
 //    double t0rms = (*tzero).second.t0rms;
 //    DTWireId wireId((*tzero).first.wheelId,
 //          (*tzero).first.stationId,
 //          (*tzero).first.sectorId,
 //          (*tzero).first.slId,
 //          (*tzero).first.layerId,
 //          (*tzero).first.cellId);
     int channelId = tzero->channelId;
     if ( channelId == 0 ) continue;
     DTWireId wireId( channelId );
     DTChamberId chamberId(wireId.chamberId());
     //double t0mean = (tzero->t0mean) - (sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
     //double t0rms = tzero->t0rms;
 // @@@ better DTT0 usage
     float t0mean_f;
     float t0rms_f;
     t0s->get(wireId,t0mean_f,t0rms_f,DTTimeUnits::counts);
     double t0mean = t0mean_f - (sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
     double t0rms = t0rms_f;
 // @@@ NEW DTT0 END
     t0sWRTChamber->set(wireId,
                t0mean,
                t0rms,
                DTTimeUnits::counts);
     if(debug){
       //cout<<"Chamber "<<chamberId<<" has reference "<<(sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
 //      cout<<"Changing t0 of wire "<<wireId<<" from "<<(*tzero).second.t0mean<<" to "<<t0mean<<endl;
       cout<<"Changing t0 of wire "<<wireId<<" from "<<tzero->t0mean<<" to "<<t0mean<<endl;
     }
   }
 
   if(debug) 
    cout << "[DTT0CalibrationNew]Writing values in DB!" << endl;
   // FIXME: to be read from cfg?
   string t0Record = "DTT0Rcd";
   // Write the t0 map to DB
   DTCalibDBUtils::writeToDB(t0Record, t0sWRTChamber);
 }

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().

                                                                  {
   string histoName;
   stringstream theStream;
   theStream << "Ch_" << wId.wheel() << "_" << wId.station() << "_" << wId.sector()
         << "_SL" << wId.superlayer() << "_L" << wId.layer() << "_W"<< wId.wire() <<"_hT0Histo";
   theStream >> histoName;
   return histoName;
 }

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().

                                                                   {
   string histoName;
   stringstream theStream;
   theStream << "Ch_" << lId.wheel() << "_" << lId.station() << "_" << lId.sector()
         << "_SL" << lId.superlayer() << "_L" << lId.layer() <<"_hT0Histo";
   theStream >> histoName;
   return histoName;
 }

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

Definition at line 54 of file DTT0CalibrationNew.h.

Referenced by util.rrapi.RRApi::dprint(), rrapi.RRApi::dprint(), pkg.AbstractPkg::generate(), util.rrapi.RRApi::get(), rrapi.RRApi::get(), pkg.AbstractPkg::get_kwds(), and pkg.AbstractPkg::write().

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.

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation