CMS 3D CMS Logo

DTT0CalibrationNew Class Reference

Analyzer class computes the mean and RMS of t0 from pulses. More...

#include <CalibMuon/DTCalibration/plugins/DTT0CalibrationNew.h>

Inheritance diagram for DTT0CalibrationNew:

edm::EDAnalyzer

List of all members.

Public Member Functions

void analyze (const edm::Event &event, const edm::EventSetup &eventSetup)
 Fill the maps with t0 (by channel).
 DTT0CalibrationNew (const edm::ParameterSet &pset)
 Constructor.
void endJob ()
 Compute the mean and the RMS of the t0 from the maps and write them to the DB with channel granularity.
virtual ~DTT0CalibrationNew ()
 Destructor.

Private Member Functions

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

Private Attributes

std::vector< std::string > cellsWithHistos
bool debug
std::string digiLabel
edm::ESHandle< DTGeometrydtGeom
unsigned int eventsForLayerT0
unsigned int eventsForWireT0
TH1D * hT0SectorHisto
std::map< DTWireId, double > mK
std::map< DTWireId, double > mK_ref
std::map< DTWireId, intnDigiPerWire
std::map< DTWireId, intnDigiPerWire_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, inttheCountT0ByChamber
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


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

Date
2008/09/01 13:29:28
Revision
1.2
Author:
S. Bolognesi - INFN Torino

Definition at line 31 of file DTT0CalibrationNew.h.


Constructor & Destructor Documentation

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

Constructor.

Definition at line 36 of file DTT0CalibrationNew.cc.

References cellsWithHistos, GenMuonPlsPt100GeV_cfg::cout, debug, digiLabel, lat::endl(), eventsForLayerT0, eventsForWireT0, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), hT0SectorHisto, nevents, rejectDigiFromPeak, retryForLayerT0, selSector, selWheel, sl, spectrum, theCalibSector, theCalibWheel, theFile, timeBoxWidth, tpPeakWidth, tpPeakWidthPerLayer, muonGeometry::wheel, and wireIdWithHistos.

00036                                                                   {
00037   // Get the debug parameter for verbose output
00038   debug = pset.getUntrackedParameter<bool>("debug");
00039   if(debug) 
00040     cout << "[DTT0CalibrationNew]Constructor called!" << endl;
00041 
00042   // Get the label to retrieve digis from the event
00043   digiLabel = pset.getUntrackedParameter<string>("digiLabel");
00044 
00045   // The root file which contain the histos per layer
00046   string rootFileName = pset.getUntrackedParameter<string>("rootFileName","DTT0PerLayer.root");
00047   theFile = new TFile(rootFileName.c_str(), "RECREATE");
00048  
00049   theCalibWheel =  pset.getUntrackedParameter<string>("calibWheel", "All"); //FIXME amke a vector of integer instead of a string
00050   if(theCalibWheel != "All") {
00051     stringstream linestr;
00052     int selWheel;
00053     linestr << theCalibWheel;
00054     linestr >> selWheel;
00055     cout << "[DTT0CalibrationNewPerLayer] chosen wheel " << selWheel << endl;
00056   }
00057 
00058   // Sector/s to calibrate
00059   theCalibSector =  pset.getUntrackedParameter<string>("calibSector", "All"); //FIXME amke a vector of integer instead of a string
00060   if(theCalibSector != "All") {
00061     stringstream linestr;
00062     int selSector;
00063     linestr << theCalibSector;
00064     linestr >> selSector;
00065     cout << "[DTT0CalibrationNewPerLayer] chosen sector " << selSector << endl;
00066   }
00067 
00068   vector<string> defaultCell;
00069   defaultCell.push_back("None");
00070   cellsWithHistos = pset.getUntrackedParameter<vector<string> >("cellsWithHisto", defaultCell);
00071   for(vector<string>::const_iterator cell = cellsWithHistos.begin(); cell != cellsWithHistos.end(); cell++){
00072     if((*cell) != "None"){
00073       stringstream linestr;
00074       int wheel,sector,station,sl,layer,wire;
00075       linestr << (*cell);
00076       linestr >> wheel >> sector >> station >> sl >> layer >> wire;
00077       wireIdWithHistos.push_back(DTWireId(wheel,station,sector,sl,layer,wire));
00078     }
00079   }
00080 
00081   hT0SectorHisto=0;
00082 
00083   nevents=0;
00084   eventsForLayerT0 = pset.getParameter<unsigned int>("eventsForLayerT0");
00085   eventsForWireT0 = pset.getParameter<unsigned int>("eventsForWireT0");
00086   tpPeakWidth = pset.getParameter<double>("tpPeakWidth");
00087   tpPeakWidthPerLayer = pset.getParameter<double>("tpPeakWidthPerLayer");
00088   timeBoxWidth = pset.getParameter<unsigned int>("timeBoxWidth"); 
00089   rejectDigiFromPeak = pset.getParameter<unsigned int>("rejectDigiFromPeak"); 
00090 
00091   spectrum = new TSpectrum(5);
00092   retryForLayerT0 = 0;  
00093 }

DTT0CalibrationNew::~DTT0CalibrationNew (  )  [virtual]

Destructor.

Definition at line 96 of file DTT0CalibrationNew.cc.

References GenMuonPlsPt100GeV_cfg::cout, debug, lat::endl(), spectrum, and theFile.

00096                                        {
00097   if(debug) 
00098     cout << "[DTT0CalibrationNew]Destructor called!" << endl;
00099 
00100   delete spectrum;
00101   theFile->Close();
00102 }


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(), DTSuperLayerId::chamberId(), GenMuonPlsPt100GeV_cfg::cout, debug, digiLabel, dtGeom, lat::endl(), edm::EventID::event(), eventsForLayerT0, eventsForWireT0, find(), edm::EventSetup::get(), getHistoName(), edm::Event::id(), it, mean(), mK, mK_ref, nDigiPerWire, nDigiPerWire_ref, nevents, qK, rejectDigiFromPeak, edm::EventID::run(), DTChamberId::sector(), selSector, selWheel, python::multivaluedict::sort(), spectrum, DTLayerId::superlayerId(), theAbsoluteT0PerWire, theCalibSector, theCalibWheel, theCountT0ByChamber, theFile, theHistoLayerMap, theHistoWireMap, theSumT0ByChamber, theTPPeakMap, timeBoxWidth, tpPeakWidth, tpPeakWidthPerLayer, DTChamberId::wheel(), and wireIdWithHistos.

00105                                                                                         {
00106   if(debug || event.id().event() % 500==0)
00107     cout << "--- [DTT0CalibrationNew] Analysing Event: #Run: " << event.id().run()
00108          << " #Event: " << event.id().event() << endl;
00109   nevents++;
00110 
00111   // Get the digis from the event
00112   Handle<DTDigiCollection> digis; 
00113   event.getByLabel(digiLabel, digis);
00114 
00115   // Get the DT Geometry
00116   eventSetup.get<MuonGeometryRecord>().get(dtGeom);
00117 
00118   // Get ttrig DB
00119   edm::ESHandle<DTTtrig> tTrigMap;
00120   eventSetup.get<DTTtrigRcd>().get(tTrigMap);
00121   
00122   // Iterate through all digi collections ordered by LayerId   
00123   DTDigiCollection::DigiRangeIterator dtLayerIt;
00124   for (dtLayerIt = digis->begin();
00125        dtLayerIt != digis->end();
00126        ++dtLayerIt){
00127     // Get the iterators over the digis associated with this LayerId
00128     const DTDigiCollection::Range& digiRange = (*dtLayerIt).second;
00129   
00130     // Get the layerId
00131     const DTLayerId layerId = (*dtLayerIt).first; //FIXME: check to be in the right sector
00132     const DTChamberId chamberId = layerId.superlayerId().chamberId();
00133 
00134     if((theCalibWheel != "All") && (layerId.superlayerId().chamberId().wheel() != selWheel))
00135       continue;
00136     if((theCalibSector != "All") && (layerId.superlayerId().chamberId().sector() != selSector))
00137       continue;
00138  
00139     //if(debug) {
00140     //  cout << "Layer " << layerId<<" with "<<distance(digiRange.first, digiRange.second)<<" digi"<<endl;
00141     //}
00142 
00143     float tTrig,tTrigRMS;
00144     tTrigMap->slTtrig(layerId.superlayerId(), tTrig, tTrigRMS);
00145     if(debug&&(nevents <= 1)){
00146         cout << "  Superlayer: " << layerId.superlayerId() << endl 
00147              << "            tTrig,tTrigRMS= " << tTrig << ", " << tTrigRMS << endl;
00148     }   
00149 
00150     // Loop over all digis in the given layer
00151     for (DTDigiCollection::const_iterator digi = digiRange.first;
00152          digi != digiRange.second;
00153          digi++) {
00154       double t0 = (*digi).countsTDC();
00155 
00156       //Use first bunch of events to fill t0 per layer
00157       if(nevents < eventsForLayerT0){
00158         //Get the per-layer histo from the map
00159         TH1I *hT0LayerHisto = theHistoLayerMap[layerId];
00160         //If it doesn't exist, book it
00161         if(hT0LayerHisto == 0){
00162           theFile->cd();
00163           float hT0Min = tTrig - 2*tTrigRMS;
00164           float hT0Max = hT0Min + timeBoxWidth;
00165           hT0LayerHisto = new TH1I(getHistoName(layerId).c_str(),
00166                                    "T0 from pulses by layer (TDC counts, 1 TDC count = 0.781 ns)",
00167                                    timeBoxWidth,hT0Min,hT0Max);
00168           if(debug)
00169             cout << "  New T0 per Layer Histo: " << hT0LayerHisto->GetName() << endl;
00170           theHistoLayerMap[layerId] = hT0LayerHisto;
00171         }
00172     
00173         //Fill the histos
00174         theFile->cd();
00175         if(hT0LayerHisto != 0) {
00176           //  if(debug)
00177           // cout<<"Filling histo "<<hT0LayerHisto->GetName()<<" with digi "<<t0<<" TDC counts"<<endl;
00178           hT0LayerHisto->Fill(t0);
00179         }
00180       }
00181 
00182       //Use all the remaining events to compute t0 per wire
00183       if(nevents>eventsForLayerT0){
00184         // Get the wireId
00185         const DTWireId wireId(layerId, (*digi).wire());
00186         if(debug) {
00187           cout << "   Wire: " << wireId << endl
00188                << "       time (TDC counts): " << (*digi).countsTDC()<< endl;
00189         }   
00190 
00191         //Fill the histos per wire for the chosen cells
00192         vector<DTWireId>::iterator it = find(wireIdWithHistos.begin(),wireIdWithHistos.end(),wireId);
00193         if (it!=wireIdWithHistos.end()){
00194           //Get the per-wire histo from the map
00195           TH1I *hT0WireHisto = theHistoWireMap[wireId]; 
00196           //If it doesn't exist, book it
00197           if(hT0WireHisto == 0){
00198             theFile->cd(); 
00199             hT0WireHisto = new TH1I(getHistoName(wireId).c_str(),"T0 from pulses by wire (TDC counts, 1 TDC count = 0.781 ns)",7000,0,7000);
00200             //hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())-100,
00201             //hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())+100);
00202             if(debug)
00203               cout << "  New T0 per wire Histo: " << hT0WireHisto->GetName() << endl;
00204             theHistoWireMap[wireId] = hT0WireHisto;
00205           }
00206           //Fill the histos
00207           theFile->cd();
00208           if(hT0WireHisto != 0) {
00209             //if(debug)
00210             // cout<<"Filling histo "<<hT0WireHisto->GetName()<<" with digi "<<t0<<" TDC counts"<<endl;
00211             hT0WireHisto->Fill(t0);
00212           }
00213         }
00214 
00215         //Check the tzero has reasonable value
00216         //float hT0Min = tTrig - 2*tTrigRMS;
00217         //float hT0Max = hT0Min + timeBoxWidth;
00218         /*if(abs(hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin()) - t0) > rejectDigiFromPeak){
00219           if(debug)
00220             cout<<"digi skipped because t0 per sector "<<hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())<<endl;
00221           continue;
00222         }*/
00223         /*if((t0 < hT0Min)||(t0 > hT0Max)){
00224           if(debug)
00225             cout<<"digi skipped because t0 outside of interval (" << hT0Min << "," << hT0Max << ")" <<endl;
00226           continue;
00227         }*/
00228         //Select per layer
00229         if(fabs(theTPPeakMap[layerId] - t0) > rejectDigiFromPeak){
00230           if(debug)
00231             cout<<"digi skipped because t0 too far from peak " << theTPPeakMap[layerId] << endl;
00232           continue;     
00233         }
00234 
00235         //Find to ref. per chamber
00236         theSumT0ByChamber[chamberId] = theSumT0ByChamber[chamberId] + t0;
00237         theCountT0ByChamber[chamberId]++;
00238 
00239         //Use second bunch of events to compute a t0 reference per wire
00240         if(nevents< (eventsForLayerT0 + eventsForWireT0)){
00241           if(!nDigiPerWire_ref[wireId]){
00242             mK_ref[wireId] = 0;
00243           }
00244           nDigiPerWire_ref[wireId] = nDigiPerWire_ref[wireId] + 1;
00245           mK_ref[wireId] = mK_ref[wireId] + (t0-mK_ref[wireId])/nDigiPerWire_ref[wireId];
00246         }
00247         //Use last all the remaining events to compute the mean and sigma t0 per wire
00248         else if(nevents>(eventsForLayerT0 + eventsForWireT0)){
00249           if(abs(t0-mK_ref[wireId]) > tpPeakWidth)
00250             continue;
00251           if(!nDigiPerWire[wireId]){
00252             theAbsoluteT0PerWire[wireId] = 0;
00253             qK[wireId] = 0;
00254             mK[wireId] = 0;
00255           }
00256           nDigiPerWire[wireId] = nDigiPerWire[wireId] + 1;
00257           theAbsoluteT0PerWire[wireId] = theAbsoluteT0PerWire[wireId] + t0;
00258           //theSigmaT0PerWire[wireId] = theSigmaT0PerWire[wireId] + (t0*t0);
00259           qK[wireId] = qK[wireId] + ((nDigiPerWire[wireId]-1)*(t0-mK[wireId])*(t0-mK[wireId])/nDigiPerWire[wireId]);
00260           mK[wireId] = mK[wireId] + (t0-mK[wireId])/nDigiPerWire[wireId];
00261         }
00262       }//end if(nevents>1000)
00263     }//end loop on digi
00264   }//end loop on layer
00265 
00266   //Use the t0 per layer histos to have an indication about the t0 position 
00267   if(nevents == eventsForLayerT0){
00268     bool increaseEvtsForLayerT0 = false;        
00269     for(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin();
00270         lHisto != theHistoLayerMap.end();
00271         lHisto++){
00272       if(debug)
00273         cout<<"Reading histogram "<<(*lHisto).second->GetName()<<" with mean "<<(*lHisto).second->GetMean()<<" and RMS "<<(*lHisto).second->GetRMS() << endl;
00274 
00275       //Find peaks
00276       //int npeaks = spectrum->Search((*lHisto).second,0.5,"goff");
00277       //int npeaks = spectrum->Search((*lHisto).second,(tpPeakWidthPerLayer/2.),"goff",0.3);
00278       int npeaks = spectrum->Search((*lHisto).second,(tpPeakWidthPerLayer/2.),"",0.3);
00279 
00280       float *peaks = spectrum->GetPositionX();  
00281       //Put in a std::vector<float>
00282       vector<float> peakMeans(peaks,peaks + npeaks);
00283       //Sort the peaks in ascending order
00284       sort(peakMeans.begin(),peakMeans.end());
00285                                 
00286       //Find best peak -- preliminary criteria: find peak closest to center of time box 
00287       float tTrig,tTrigRMS;
00288       tTrigMap->slTtrig((*lHisto).first.superlayerId(), tTrig, tTrigRMS);
00289       float timeBoxCenter = (2*tTrig + (float)timeBoxWidth)/2.; 
00290       float hMin = (*lHisto).second->GetXaxis()->GetXmin();
00291       float hMax = (*lHisto).second->GetXaxis()->GetXmax();             
00292       vector<float>::const_iterator tpPeak = peakMeans.end();
00293       for(vector<float>::const_iterator it = peakMeans.begin(); it != peakMeans.end(); ++it){
00294         float mean = *it;
00295 
00296         int bin = (*lHisto).second->GetXaxis()->FindBin(mean);
00297         float yp = (*lHisto).second->GetBinContent(bin);
00298         if(debug) cout << "Peak : (" << mean << "," << yp << ")" << endl; 
00299 
00300         //Find RMS
00301         float previous_peak = (it == peakMeans.begin())?hMin:*(it - 1);
00302         float next_peak = (it == (peakMeans.end()-1))?hMax:*(it + 1);
00303 
00304         float rangemin = mean - (mean - previous_peak)/8.;
00305         float rangemax = mean + (next_peak - mean)/8.;
00306         int binmin = (*lHisto).second->GetXaxis()->FindBin(rangemin);
00307         int binmax = (*lHisto).second->GetXaxis()->FindBin(rangemax);
00308         (*lHisto).second->GetXaxis()->SetRange(binmin,binmax);
00309         //RMS estimate
00310         float rms_seed = (*lHisto).second->GetRMS();
00311 
00312         /*rangemin = mean - 2*rms_seed;
00313         rangemax = mean + 2*rms_seed;
00314         if(debug) cout << "Seed for RMS, Fit min, Fit max: " << rms_seed << ", " << rangemin << ", " << rangemax << endl;
00315         //Fit to gaussian
00316         string funcname("fitFcn_");
00317         funcname += (*lHisto).second->GetName();
00318         if(debug) cout << "Fitting function " << funcname << endl; 
00319         TF1* func = new TF1(funcname.c_str(),"gaus",rangemin,rangemax);
00320         func->SetParameters(yp,mean,rms_seed);
00321         (*lHisto).second->Fit(func,"Q","",rangemin,rangemax);
00322 
00323         float fitconst = func->GetParameter(0);
00324         float fitmean = func->GetParameter(1);
00325         float fitrms = func->GetParameter(2);
00326         float chisquare = func->GetChisquare()/func->GetNDF();
00327         if(debug) cout << "Gaussian fit constant,mean,RMS,chi2= " << fitconst << ", " << fitmean << ", " << fitrms << ", " << chisquare << endl;*/
00328 
00329         //Reject peaks with RMS larger than specified
00330         //if(fitrms > tpPeakWidth) continue;
00331         if(rms_seed > tpPeakWidthPerLayer) continue;
00332 
00333         if(fabs(mean - timeBoxCenter) < fabs(*tpPeak - timeBoxCenter)) tpPeak = it;
00334       } 
00335       //Didn't find peak        
00336       /*if(tpPeak == peakMeans.end()){
00337         if(retryForLayerT0 < 2){
00338           increaseEvtsForLayerT0 = true;
00339           retryForLayerT0++;
00340           break;
00341         } 
00342       }*/
00343 
00344       float selPeak = (tpPeak != peakMeans.end())?*tpPeak:(*lHisto).second->GetBinCenter((*lHisto).second->GetMaximumBin());            
00345       if(debug) cout << "Peak selected at " << selPeak << endl;
00346         
00347       theTPPeakMap[(*lHisto).first] = selPeak;
00348                 
00349       //Take the mean as a first t0 estimation
00350       /*if((*lHisto).second->GetRMS() < tpPeakWidth){
00351         if(hT0SectorHisto == 0){
00352           hT0SectorHisto = new TH1D("hT0AllLayerOfSector","T0 from pulses per layer in sector", 
00353                                     //20, (*lHisto).second->GetMean()-100, (*lHisto).second->GetMean()+100);
00354                                     700, 0, 7000);
00355                                     //300,3300,3600);   
00356         }
00357         if(debug)
00358           cout<<" accepted"<<endl;
00359         //TH1I* aux_histo = (*lHisto).second;
00360         //aux_histo->GetXaxis()->SetRangeUser(3300,3600);
00361         hT0SectorHisto->Fill((*lHisto).second->GetMean());
00362         //hT0SectorHisto->Fill(aux_histo->GetMean());
00363       }
00364       //Take the mean of noise + 400ns as a first t0 estimation
00365       //if((*lHisto).second->GetRMS()>10.0 && ((*lHisto).second->GetRMS()<15.0)){
00366       //double t0_estim = (*lHisto).second->GetMean() + 400;
00367       //if(hT0SectorHisto == 0){
00368       //  hT0SectorHisto = new TH1D("hT0AllLayerOfSector","T0 from pulses per layer in sector", 
00369       //                            //20, t0_estim-100, t0_estim+100);
00370       //                            700, 0, 7000);
00371       //}
00372       //if(debug)
00373       //  cout<<" accepted + 400ns"<<endl;
00374       //hT0SectorHisto->Fill((*lHisto).second->GetMean() + 400);
00375       //}
00376       if(debug)
00377         cout<<endl;
00378 
00379       theT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetMean();
00380       theSigmaT0LayerMap[(*lHisto).second->GetName()] = (*lHisto).second->GetRMS();*/
00381     }
00382     /*if(!hT0SectorHisto){
00383       cout<<"[DTT0CalibrationNew]: All the t0 per layer are still uncorrect: trying with greater number of events"<<endl;
00384       eventsForLayerT0 = eventsForLayerT0*2;
00385       return;
00386     }
00387     if(debug)
00388       cout<<"[DTT0CalibrationNew] t0 reference for this sector "<<
00389         hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin())<<endl;*/
00390     if(increaseEvtsForLayerT0){
00391         cout<<"[DTT0CalibrationNew]: t0 per layer are still uncorrect: trying with greater number of events"<<endl;
00392         eventsForLayerT0 = eventsForLayerT0*2;
00393         return;
00394     }           
00395   } 
00396 }

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

References funct::abs(), DTT0::begin(), DTTimeUnits::counts, GenMuonPlsPt100GeV_cfg::cout, debug, dtGeom, DTT0::end(), lat::endl(), nDigiPerWire, qK, DTT0::set(), sl, funct::sqrt(), theAbsoluteT0PerWire, theCountT0ByChamber, theFile, theHistoLayerMap, theHistoWireMap, theRefT0ByChamber, theRelativeT0PerWire, theSigmaT0PerWire, theSumT0ByChamber, tzero, and DTCalibDBUtils::writeToDB().

00398                                 {
00399 
00400   DTT0* t0s = new DTT0();
00401   DTT0* t0sWRTChamber = new DTT0();
00402 
00403   if(debug) 
00404     cout << "[DTT0CalibrationNewPerLayer]Writing histos to file!" << endl;
00405 
00406   theFile->cd();
00407   //hT0SectorHisto->Write();
00408   for(map<DTWireId, TH1I*>::const_iterator wHisto = theHistoWireMap.begin();
00409       wHisto != theHistoWireMap.end();
00410       wHisto++) {
00411     (*wHisto).second->Write(); 
00412   }
00413   for(map<DTLayerId, TH1I*>::const_iterator lHisto = theHistoLayerMap.begin();
00414       lHisto != theHistoLayerMap.end();
00415       lHisto++) {
00416     (*lHisto).second->Write(); 
00417   }  
00418 
00419   if(debug) 
00420     cout << "[DTT0CalibrationNew] Compute and store t0 and sigma per wire" << endl;
00421 
00422   for(map<DTChamberId,double>::const_iterator chamber = theSumT0ByChamber.begin();
00423       chamber != theSumT0ByChamber.end();
00424       ++chamber) theRefT0ByChamber[(*chamber).first] = theSumT0ByChamber[(*chamber).first]/((double)theCountT0ByChamber[(*chamber).first]);
00425   
00426   for(map<DTWireId, double>::const_iterator wiret0 = theAbsoluteT0PerWire.begin();
00427       wiret0 != theAbsoluteT0PerWire.end();
00428       wiret0++){
00429     if(nDigiPerWire[(*wiret0).first]){
00430       double t0 = (*wiret0).second/nDigiPerWire[(*wiret0).first];
00431       DTChamberId chamberId = ((*wiret0).first).chamberId();
00432       //theRelativeT0PerWire[(*wiret0).first] = t0 - hT0SectorHisto->GetBinCenter(hT0SectorHisto->GetMaximumBin());
00433       theRelativeT0PerWire[(*wiret0).first] = t0 - theRefT0ByChamber[chamberId];        
00434       cout<<"Wire "<<(*wiret0).first<<" has    t0 "<<t0<<"(absolute) "<<theRelativeT0PerWire[(*wiret0).first]<<"(relative)";
00435 
00436       //theSigmaT0PerWire[(*wiret0).first] = sqrt((theSigmaT0PerWire[(*wiret0).first] / nDigiPerWire[(*wiret0).first]) - t0*t0);
00437       theSigmaT0PerWire[(*wiret0).first] = sqrt(qK[(*wiret0).first]/nDigiPerWire[(*wiret0).first]);
00438       cout<<"    sigma "<<theSigmaT0PerWire[(*wiret0).first]<<endl;
00439     }
00440     else{
00441       cout<<"[DTT0CalibrationNew] ERROR: no digis in wire "<<(*wiret0).first<<endl;
00442       abort();
00443     }
00444   }
00445 
00447   // Get all the sls from the setup
00448   const vector<DTSuperLayer*> superLayers = dtGeom->superLayers();     
00449   // Loop over all SLs
00450   for(vector<DTSuperLayer*>::const_iterator  sl = superLayers.begin();
00451       sl != superLayers.end(); sl++) {
00452 
00453   
00454     //Compute mean for odd and even superlayers
00455     double oddLayersMean=0;
00456     double evenLayersMean=0; 
00457     double oddLayersDen=0;
00458     double evenLayersDen=0;
00459     for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
00460         wiret0 != theRelativeT0PerWire.end();
00461         wiret0++){
00462       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
00463         if(debug)
00464           cout<<"[DTT0CalibrationNew] Superlayer "<<(*sl)->id()
00465               <<"layer " <<(*wiret0).first.layerId().layer()<<" with "<<(*wiret0).second<<endl;
00466         if(((*wiret0).first.layerId().layer()) % 2){
00467           oddLayersMean = oddLayersMean + (*wiret0).second;
00468           oddLayersDen++;
00469         }
00470         else{
00471           evenLayersMean = evenLayersMean + (*wiret0).second;
00472           evenLayersDen++;
00473         }
00474       }
00475     }
00476     oddLayersMean = oddLayersMean/oddLayersDen;
00477     evenLayersMean = evenLayersMean/evenLayersDen;
00478     if(debug && oddLayersMean)
00479       cout<<"[DTT0CalibrationNew] Relative T0 mean for  odd layers "<<oddLayersMean<<"  even layers"<<evenLayersMean<<endl;
00480 
00481     //Compute sigma for odd and even superlayers
00482     double oddLayersSigma=0;
00483     double evenLayersSigma=0;
00484     for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
00485         wiret0 != theRelativeT0PerWire.end();
00486         wiret0++){
00487       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
00488         if(((*wiret0).first.layerId().layer()) % 2){
00489           oddLayersSigma = oddLayersSigma + ((*wiret0).second - oddLayersMean) * ((*wiret0).second - oddLayersMean);
00490         }
00491         else{
00492           evenLayersSigma = evenLayersSigma + ((*wiret0).second - evenLayersMean) * ((*wiret0).second - evenLayersMean);
00493         }
00494       }
00495     }
00496     oddLayersSigma = oddLayersSigma/oddLayersDen;
00497     evenLayersSigma = evenLayersSigma/evenLayersDen;
00498     oddLayersSigma = sqrt(oddLayersSigma);
00499     evenLayersSigma = sqrt(evenLayersSigma);
00500 
00501     if(debug && oddLayersMean)
00502       cout<<"[DTT0CalibrationNew] Relative T0 sigma for  odd layers "<<oddLayersSigma<<"  even layers"<<evenLayersSigma<<endl;
00503 
00504     //Recompute the mean for odd and even superlayers discarding fluctations
00505     double oddLayersFinalMean=0; 
00506     double evenLayersFinalMean=0;
00507     for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
00508         wiret0 != theRelativeT0PerWire.end();
00509         wiret0++){
00510       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
00511         if(((*wiret0).first.layerId().layer()) % 2){
00512           if(abs((*wiret0).second - oddLayersMean) < (2*oddLayersSigma))
00513             oddLayersFinalMean = oddLayersFinalMean + (*wiret0).second;
00514         }
00515         else{
00516           if(abs((*wiret0).second - evenLayersMean) < (2*evenLayersSigma))
00517             evenLayersFinalMean = evenLayersFinalMean + (*wiret0).second;
00518         }
00519       }
00520     }
00521     oddLayersFinalMean = oddLayersFinalMean/oddLayersDen;
00522     evenLayersFinalMean = evenLayersFinalMean/evenLayersDen;
00523     if(debug && oddLayersMean)
00524       cout<<"[DTT0CalibrationNew] Final relative T0 mean for  odd layers "<<oddLayersFinalMean<<"  even layers"<<evenLayersFinalMean<<endl;
00525 
00526     for(map<DTWireId, double>::const_iterator wiret0 = theRelativeT0PerWire.begin();
00527         wiret0 != theRelativeT0PerWire.end();
00528         wiret0++){
00529       if((*wiret0).first.layerId().superlayerId() == (*sl)->id()){
00530         double t0=-999;
00531         if(((*wiret0).first.layerId().layer()) % 2)
00532           t0 = (*wiret0).second + (evenLayersFinalMean - oddLayersFinalMean);
00533         else
00534           t0 = (*wiret0).second;
00535 
00536         cout<<"[DTT0CalibrationNew] Wire "<<(*wiret0).first<<" has    t0 "<<(*wiret0).second<<" (relative, after even-odd layer corrections)  "
00537             <<"    sigma "<<theSigmaT0PerWire[(*wiret0).first]<<endl;
00538         //Store the results into DB
00539         t0s->set((*wiret0).first, t0, theSigmaT0PerWire[(*wiret0).first],DTTimeUnits::counts); 
00540       }
00541     }
00542   }
00543   
00545   if(debug) 
00546     cout << "[DTT0CalibrationNew]Computing relative t0 wrt to chamber average" << endl;
00547   //Compute the reference for each chamber
00548   map<DTChamberId,double> sumT0ByChamber;
00549   map<DTChamberId,int> countT0ByChamber;
00550   for(DTT0::const_iterator tzero = t0s->begin();
00551       tzero != t0s->end(); tzero++) {
00552     DTChamberId chamberId((*tzero).first.wheelId,
00553                           (*tzero).first.stationId,
00554                           (*tzero).first.sectorId);
00555     sumT0ByChamber[chamberId] = sumT0ByChamber[chamberId] + (*tzero).second.t0mean;
00556     countT0ByChamber[chamberId]++;
00557   }
00558 
00559   //Change reference for each wire and store the new t0s in the new map
00560   for(DTT0::const_iterator tzero = t0s->begin();
00561       tzero != t0s->end(); tzero++) {
00562     DTChamberId chamberId((*tzero).first.wheelId,
00563                           (*tzero).first.stationId,
00564                           (*tzero).first.sectorId);
00565     double t0mean = ((*tzero).second.t0mean) - (sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
00566     double t0rms = (*tzero).second.t0rms;
00567     DTWireId wireId((*tzero).first.wheelId,
00568                     (*tzero).first.stationId,
00569                     (*tzero).first.sectorId,
00570                     (*tzero).first.slId,
00571                     (*tzero).first.layerId,
00572                     (*tzero).first.cellId);
00573     t0sWRTChamber->set(wireId,
00574                        t0mean,
00575                        t0rms,
00576                        DTTimeUnits::counts);
00577     if(debug){
00578       //cout<<"Chamber "<<chamberId<<" has reference "<<(sumT0ByChamber[chamberId]/countT0ByChamber[chamberId]);
00579       cout<<"Changing t0 of wire "<<wireId<<" from "<<(*tzero).second.t0mean<<" to "<<t0mean<<endl;
00580     }
00581   }
00582 
00584   if(debug) 
00585    cout << "[DTT0CalibrationNew]Writing values in DB!" << endl;
00586   // FIXME: to be read from cfg?
00587   string t0Record = "DTT0Rcd";
00588   // Write the t0 map to DB
00589   DTCalibDBUtils::writeToDB(t0Record, t0sWRTChamber);
00590 }

string DTT0CalibrationNew::getHistoName ( const DTLayerId lId  )  const [private]

Definition at line 601 of file DTT0CalibrationNew.cc.

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

00601                                                                   {
00602   string histoName;
00603   stringstream theStream;
00604   theStream << "Ch_" << lId.wheel() << "_" << lId.station() << "_" << lId.sector()
00605             << "_SL" << lId.superlayer() << "_L" << lId.layer() <<"_hT0Histo";
00606   theStream >> histoName;
00607   return histoName;
00608 }

string DTT0CalibrationNew::getHistoName ( const DTWireId wId  )  const [private]

Definition at line 592 of file DTT0CalibrationNew.cc.

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

Referenced by analyze().

00592                                                                  {
00593   string histoName;
00594   stringstream theStream;
00595   theStream << "Ch_" << wId.wheel() << "_" << wId.station() << "_" << wId.sector()
00596             << "_SL" << wId.superlayer() << "_L" << wId.layer() << "_W"<< wId.wire() <<"_hT0Histo";
00597   theStream >> histoName;
00598   return histoName;
00599 }


Member Data Documentation

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

Definition at line 102 of file DTT0CalibrationNew.h.

Referenced by DTT0CalibrationNew().

bool DTT0CalibrationNew::debug [private]

Definition at line 56 of file DTT0CalibrationNew.h.

Referenced by analyze(), DTT0CalibrationNew(), endJob(), and ~DTT0CalibrationNew().

std::string DTT0CalibrationNew::digiLabel [private]

Definition at line 59 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

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

Definition at line 125 of file DTT0CalibrationNew.h.

Referenced by analyze(), and endJob().

unsigned int DTT0CalibrationNew::eventsForLayerT0 [private]

Definition at line 69 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

unsigned int DTT0CalibrationNew::eventsForWireT0 [private]

Definition at line 71 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

TH1D* DTT0CalibrationNew::hT0SectorHisto [private]

Definition at line 96 of file DTT0CalibrationNew.h.

Referenced by DTT0CalibrationNew().

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

Definition at line 110 of file DTT0CalibrationNew.h.

Referenced by analyze().

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

Definition at line 111 of file DTT0CalibrationNew.h.

Referenced by analyze().

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

Definition at line 108 of file DTT0CalibrationNew.h.

Referenced by analyze(), and endJob().

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

Definition at line 109 of file DTT0CalibrationNew.h.

Referenced by analyze().

unsigned int DTT0CalibrationNew::nevents [private]

Definition at line 67 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

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

Definition at line 112 of file DTT0CalibrationNew.h.

Referenced by analyze(), and endJob().

unsigned int DTT0CalibrationNew::rejectDigiFromPeak [private]

Definition at line 83 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

unsigned int DTT0CalibrationNew::retryForLayerT0 [private]

Definition at line 85 of file DTT0CalibrationNew.h.

Referenced by DTT0CalibrationNew().

int DTT0CalibrationNew::selSector [private]

Definition at line 91 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

int DTT0CalibrationNew::selWheel [private]

Definition at line 89 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

TSpectrum* DTT0CalibrationNew::spectrum [private]

Definition at line 98 of file DTT0CalibrationNew.h.

Referenced by analyze(), DTT0CalibrationNew(), and ~DTT0CalibrationNew().

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

Definition at line 105 of file DTT0CalibrationNew.h.

Referenced by analyze(), and endJob().

std::string DTT0CalibrationNew::theCalibSector [private]

Definition at line 90 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

std::string DTT0CalibrationNew::theCalibWheel [private]

Definition at line 88 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

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

Definition at line 121 of file DTT0CalibrationNew.h.

Referenced by analyze(), and endJob().

TFile* DTT0CalibrationNew::theFile [private]

Definition at line 62 of file DTT0CalibrationNew.h.

Referenced by analyze(), DTT0CalibrationNew(), endJob(), and ~DTT0CalibrationNew().

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

Definition at line 94 of file DTT0CalibrationNew.h.

Referenced by analyze(), and endJob().

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

Definition at line 114 of file DTT0CalibrationNew.h.

Referenced by analyze(), and endJob().

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.

Referenced by endJob().

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

Definition at line 106 of file DTT0CalibrationNew.h.

Referenced by endJob().

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.

Referenced by endJob().

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

Definition at line 120 of file DTT0CalibrationNew.h.

Referenced by analyze(), and endJob().

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.

Referenced by analyze().

unsigned int DTT0CalibrationNew::timeBoxWidth [private]

Definition at line 80 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

double DTT0CalibrationNew::tpPeakWidth [private]

Definition at line 74 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

double DTT0CalibrationNew::tpPeakWidthPerLayer [private]

Definition at line 77 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().

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

Definition at line 101 of file DTT0CalibrationNew.h.

Referenced by analyze(), and DTT0CalibrationNew().


The documentation for this class was generated from the following files:
Generated on Tue Jun 9 18:19:07 2009 for CMSSW by  doxygen 1.5.4